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ON  THE  COMBINED  PERFORMANCE  OF  NON-LOCAL  ARTIFICIAL  BOUNDARY 
CONDITIONS  WITH  THE  NEW  GENERATION  OF  ADVANCED  MULTIGRID  FLOW 

SOLVERS  * 

T.  W.  ROBERTS*,  D.  SIDILKOVER* ,  AND  S.  V.  TSYNKOV§ 

Abstract.  We  develop  theoretically  and  implement  numerically  a  unified  flow  solution  methodology 
that  combines  the  advantages  relevant  to  two  independent  groups  of  methods  in  CFD  that  have  recently 
proven  successful:  The  new  factorizable  schemes  for  the  equations  of  hydrodynamics  that  facilitate  the 
construction  of  optimally  convergent  multigrid  algorithms,  and  highly  accurate  global  far-field  artificial 
boundary  conditions  (ABCs).  The  primary  result  that  we  have  obtained  is  the  following.  Global  ABCs  do 
not  hamper  the  optimal  (i.e.,  unimprovable)  multigrid  convergence  rate  pertinent  to  the  solver.  At  the  same 
time,  contrary  to  the  standard  local  ABCs,  the  solution  accuracy  provided  by  the  global  ABCs  deteriorates 
very  slightly  or  does  not  deteriorate  at  all  when  the  computational  domain  shrinks,  which  clearly  translates 
into  substantial  savings  of  computer  resources. 
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convergence  rate,  exact  solution,  error  profiles 
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1*  In^r°duction.  The  subject  of  this  paper  is  development  of  the  joint  formulation,  combined  im- 
plementation,  and  subsequent  performance  assessment  for  the  exact  nonlocal  far-field  artificial  boundary 
conditions  (ABCs)  coupled  with  the  new  generation  of  multigrid  flow  solvers  based  on  factorizable  schemes 
for  the  equations  of  hydrodynamics.  Both  methodologies  have  independently  proven  efficient  and  promising 
for  the  numerical  solution  of  different  flow  problems  and  thus  it  appears  natural  to  try  and  analyze  their 
performance  if  combined  with  one  another.  Accordingly,  we  organize  the  material  as  follows.  First,  we 
briefly  review  relevant  results  in  the  two  aforementioned  independent  areas:  Flow  solvers  and  ABCs,  and 
then  formulate  the  motivation  and  objectives  for  the  current  study.  Next  follows  the  core  of  the  paper  — 
description  of  the  algorithm  and  the  set  of  numerical  experiments.  Finally,  we  discuss  the  results  that  we 
have  obtained,  as  well  as  possible  extensions  for  the  future. 

1.1.  Background  on  the  scheme.  To  emphasize  the  advantages  of  the  new  generation  of  flow  solvers 
that  emerge  nowadays,  we  first  briefly  review  the  history  of  developing  the  numerical  methods  for  compress¬ 
ible  flow  computations. 
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The  Full-Potential  equation  is  one  of  the  simplest  mathematical  models  describing  the  compressible  fluid 
flow  phenomena.  Active  research  on  numerical  methods  for  solving  this  equation  has  begun  in  the  early  70’s 
with  the  work  of  Murman  and  Cole  [1].  Numerical  schemes  were  developed  for  two-  and  three-dimensional 
cases.  Some  research  has  also  been  done  on  applying  the  multigrid  methods  as  a  means  to  accelerate  the 
convergence  to  steady  state. 

The  Euler  equations  constitute  a  more  general  model  for  the  fluid  flow  since  they  account  for  the  ad- 
vection  of  vorticity  and  entropy  as  well.  The  first  rather  efficient  method  for  the  steady-state  computations 
using  Euler’s  equations  was  developed  in  early  80’s  by  Jameson  and  co-authors  [2,3];  it  still  remains  (with 
some  modifications)  widely  used  in  modern  computational  practice.  However,  there  was  no  obvious  con¬ 
nection  between  this  new  Euler  solver  and  the  previously  used  Full-Potential  solvers.  The  research  on  the 
Full-Potential  solvers  was  effectively  abandoned  in  its  infancy  in  favor  of  the  solvers  for  the  Euler  equations. 

On  the  other  hand,  it  is  becoming  increasingly  clear  nowadays  that  the  methods  currently  used  for 
the  steady  compressible  Euler  (and  Navier- Stokes)  computations  have  severe  limitations  in  terms  of  both 
computational  efficacy  and  robustness. 

A  major  difficulty  for  the  numerical  treatment  of  compressible  flow  is  the  possible  presence  of  shocks  in 
the  solution.  It  is  well  known  that  a  scheme,  which  is  both  second  order  accurate  and  avoids  undershoots 
and  overshoots  near  discontinuities  (which  may  trigger  nonlinear  instability),  has  to  be  nonlinear.  Such 
a  scheme  has  to  incorporate  the  so-called  high-resolution  mechanism,  i.e.,  a  smoothness  monitor,  that  is 
usually  implemented  in  the  form  of  a  flux  limiter.  Initially,  the  schemes  of  this  type  were  developed  for 
one-dimensional  case.  Then,  they  were  extended  to  multiple  dimensions  using  a  dimension-by-dimension  ap¬ 
proach.  However,  these  straightforward  multidimensional  extensions  of  high-resolution  schemes  have  proven 
ineffective  when  applied  for  calculating  the  steady-state  aerodynamic  solutions,  which  are  of  great  interest 
for  applications.  It  turns  out  that  the  standard  multidimensional  high-resolution  discretizations  (obtained 
by  a  dimension-by-dimension  extension)  suffer  from  the  following  deficiency:  The  contribution  of  the  high- 
frequency  error  components  to  the  residuals  of  the  discrete  equations  is  very  small,  which  translates  into 
the  poor  stability  characteristics  of  the  steady-state  discretization.  This  makes  the  construction  of  a  good 
relaxation  scheme,  in  particular,  smoother  for  a  multigrid  procedure,  inherently  difficult,  if  possible  at  all. 
For  example,  Spekreijse  [4]  had  observed  that  the  Gauss-Seidel  relaxation  was  unstable  when  implemented 
along  with  such  schemes.  Thus,  the  multigrid  methods  based  on  the  Gauss-Seidel  relaxation  could  not  be 
used  to  accelerate  the  convergence  to  steady  state.  Instead,  the  multigrid  solvers  designed  for  routine  use  had 
to  resort  to  a  defect-correction  technique  or  multistage  Runge-Kutta  relaxation,  and  as  such  their  efficiency 
was  relatively  poor. 

A  genuinely  multidimensional  advection  scheme  was  constructed  in  [5,6].  The  scheme  was  named  “gen¬ 
uinely  multidimensional”  since  it  imitates  well  the  anisotropy  of  the  advection  phenomena  in  two  dimensions: 
The  artificial  dissipation  is  added  only  along  the  streamline,  while  the  high-resolution  mechanism  affects  sig¬ 
nificantly  the  cross-flow  direction  only.  A  key  feature  of  this  scheme  is  its  two-dimensional  limiter,  i.e.,  the 
limiter- function  of  an  argument  that  is  the  ratio  of  divided  differences  in  two  different  coordinate  directions. 
The  scheme  was  formulated  in  the  control-volume  context  for  Cartesian  grids  and  relied  on  the  compact 
9-point-box  stencil.  The  fundamental  advantage  of  this  approach  is  that  the  two-dimensional  high-resolution 
mechanism  does  not  damage  the  stability  properties  of  the  steady-state  discretization. 

Generalization  of  these  ideas  to  the  systems  of  equations  was  not  straightforward  and  took  substantial 
time  and  effort.  A  generalization  of  this  type  that  yields  a  robust  scheme  for  the  Euler  equations,  which  is 
suitable  for  the  computations  of  a  wide  range  of  flow  regimes,  was  presented  in  [7].  Later  it  was  described  in 
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more  detail,  including  the  implications  for  multigrid  and  three-dimensional  extensions,  see  [8]  and  [9].  The 
key  idea  was  not  to  try  and  directly  apply  the  scalar  advection  scheme  to  the  case  of  systems,  but  rather 
start  from  the  scratch  and  use  for  the  systems  the  same  strategy  that  was  employed  to  construct  the  scalar 
scheme.  The  resulting  scheme  for  the  Euler  equations  was  shown  to  generate  a  very  good  quality  solutions 
for  subsonic,  transonic,  and  supersonic  regimes.  A  manifestation  of  “the  genuine  multidimensionality”  of 
this  scheme  was  rotationally-invariant  form  of  its  artificial  dissipation. 

As  in  the  scalar  case,  the  fundamental  advantage  of  the  approach  [7-9]  is  that  its  high-resolution  mecha¬ 
nism  does  not  damage  the  steady-state  stability  properties  of  the  scheme.  Similarly  to  the  scalar  case,  it  was 
demonstrated  in  [8,9]  that  the  Gauss-Seidel  relaxation  is  stable  when  applied  directly  to  the  resulting  high- 
resolution  discretization  of  the  hyperbolic  systems.  This  yields  a  very  simple,  efficient  and  robust  multigrid 
solver  for  the  compressible  Euler  equations  suitable  for  the  entire  range  of  flow  regimes. 

Besides  the  aforementioned  difficulties  in  the  iterative  solution  of  the  steady-state  flow  equations  that 
originate  from  the  poor  sensitivity  of  the  residuals  to  the  high-frequency  error  content,  and  that  have  been 
largely  overcome  by  introducing  the  genuinely  multidimensional  schemes,  multigrid  iterative  solvers  are 
prone  to  deficiencies  of  a  different  nature  that  can  hamper  the  performance  even  if  a  good  smoother  is 
available.  As  outlined  in  [10],  for  the  advection  dominated  problems  the  coarse  grid  provides  only  a  fraction 
of  the  needed  correction  for  certain  error  components.  Unlike  the  previous  frequency-based  description, 
this  time  the  “problematic”  contributions  to  the  error  are  the  so-called  characteristic  components,  see  [10] 
(frequency- wise,  they  are  typically  mid-range).  On  the  other  hand,  it  is  well  known  that  the  steady  Euler 
equations  can  be  factored  into  the  advection  and  Full-Potential  parts;  the  latter  is  of  either  elliptic  or 
hyperbolic  type  depending  on  the  flow  regime  (subsonic  or  supersonic,  respectively).  The  aforementioned 
difficulty  (insufficient  coarse-grid  corrections  for  characteristic  error  components)  can  be  avoided  (see  [10]) 
by  constructing  a  solver  that  distinguishes  between  the  different  factors  of  the  system  and  treats  each  one 
appropriately.  In  the  subsonic  case,  for  instance,  the  advection  factor  can  be  treated  by  marching  and  the 
elliptic  factor  —  by  multigrid.  The  efficiency  of  such  an  algorithm  will  be  essentially  the  same  as  of  a  multigrid 
solver  applied  to  the  elliptic  part  only.  Such  algorithms  are  referred  to  as  “essentially  optimal.”  An  approach 
to  separating  the  co-factors  —  the  so-called  Distributive  Gauss-Seidel  relaxation  —  was  proposed  in  [10].  It 
was  demonstrated  in  [11]  that  using  this  approach  one  can  obtain  an  essentially  optimal  multigrid  efficiency 
for  a  staggered-grid  discretization  of  the  incompressible  Navier-Stokes  equations;  a  similar  observation  was 
made  earlier  in  [12]. 

Another  way  (alternative  to  the  Distributive  Gauss-Seidel  relaxation)  toward  exploiting  the  ellip¬ 
tic/hyperbolic  distinction  in  the  governing  equations  and  thus  achieving  the  optimal  multigrid  efficiency 
is  based  on  the  well-known  pressure  formulation  of  the  Euler  equations,  which  also  amounts  to  the  factor¬ 
ization  of  the  system  into  the  elliptic  and  advection  parts.  A  particular  formulation  of  the  scheme  and  the 
corresponding  multigrid  algorithm  that  employ  this  idea  was  proposed  in  [13].  This  approach  was  further 
generalized  in  [14];  work  [14]  also  contains  an  extensive  set  of  numerical  computations  and  more  detail  re¬ 
garding  the  implementation  of  the  scheme.  The  main  advantage  of  this  approach  is  its  simplicity,  it  can 
also  be  classified  as  Weighted  Gauss-Seidel  relaxation  [10].  Subsequent  work  in  this  direction  is  presented 
in  [15,16].  The  limitation  of  this  approach,  however,  is  that  it  is  not  clear  as  of  yet  whether  or  not  it  can  be 
generalized  to  the  case  of  viscous  compressible  flows. 

It,  however,  turns  out,  see  [17],  that  the  idea  of  factorization  that  has  successfully  led  to  the  construction 
of  essentially  optimal  multigrid  solvers  in  [10-16]  cannot  be  applied  to  an  important  class  of  schemes  that  is 
routinely  used  for  compressible  flow  computations  —  the  shock-capturing  schemes  of  all  types,  including  the 
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original  formulation  of  genuinely  multidimensional  scheme,  see  [7-9].  As  these  schemes  possess  a  collection  of 
features  that  cannot  be  either  compromised  or  traded  for  in  computational  practice,  primarily  the  mechanism 
for  handling  discontinuous  solutions,  there  is  a  need  for  discretization  that  would  rather  add  the  capability  of 
factoring  the  governing  PDEs  to  the  existing  properties  of  the  scheme.  As  has  been  pointed  out  in  [17],  there 
is  a  certain  freedom  in  the  way  how  the  multidimensional  corrections  are  built  for  genuinely  multidimensional 
schemes  [7,8].  These  extra  freedom  can  be  used  to  make  the  resulting  scheme  also  factorizable ,  i.e.,  make 
it  capable  of  reflecting  the  mixed  nature  of  the  governing  PDEs.  This  means  that  the  co-factors  of  the 
different  type  can  be  distinguished  directly  at  the  discrete  level  (while  still  keeping  the  multidimensional 
shock-capturing  property).  This,  in  turn,  facilitates  the  construction  of  an  optimally  efficient  multigrid 
solver  through  the  design  of  a  special  Distributive  relaxation.  The  scheme  of  this  type  was  described  in  [18] 
with  the  emphasis  on  the  subsonic  case.  The  construction  was  later  extended  to  the  case  of  three  spatial 
dimensions  [19]  and  generalized  coordinates  [20].  Extending  the  approach  to  the  transonic  and  supersonic 
regimes  is  underway. 

1.2.  Background  on  the  ABCs.  Artificial  boundary  conditions  furnish  a  widely  used  approach  for 
the  numerical  treatment  of  boundary-value  problems  initially  formulated  on  unbounded  domains.  These 
boundary  conditions  are  typically  set  at  the  external  boundary  of  a  finite  computational  domain  once  the 
latter  is  obtained  from  the  original  unbounded  domain  by  means  of  truncation.  Implementation  of  the 
ABCs  completes  the  “truncated  problem”  and  therefore,  makes  it  available  for  solution  on  the  computer. 
Different  authors  have  repeatedly  shown  both  theoretically  and  experimentally  that  the  overall  accuracy  and 
performance  of  numerical  algorithms,  as  well  as  interpretation  of  the  results,  strongly  depend  on  the  proper 
treatment  of  external  artificial  boundaries. 

The  choice  of  the  ABCs  is  typically  not  unique.  Clearly,  the  minimal  necessary  requirement  of  the  ABCs 
is  to  ensure  the  solvability  of  the  truncated  problem.  However,  meeting  this  requirement  only  does  not  guar¬ 
antee  that  the  solution  found  inside  the  computational  domain  will  be  anywhere  close  to  the  corresponding 
fragment  of  the  solution  to  the  original  (infinite-domain)  problem.  Therefore,  we  must  additionally  require 
that  the  two  solutions  be  in  a  certain  sense  close  to  one  another  on  the  truncated  domain.  Ideally,  these  two 
solutions  coincide,  which  corresponds  to  the  so-called  exact  ABCs. 

It  turns  out  that  in  most  cases,  the  exact  ABCs  are  nonlocal ,  for  steady-state  problems  in  space  and 
for  time-dependent  problems  also  in  time.  Besides,  many  methodologies  for  obtaining  exact  ABCs  lack  geo¬ 
metric  universality  as  they  rely  on  integral  transforms  along  the  boundary  and  the  separation  of  variables. 
In  practical  computing,  the  aforementioned  nonlocality  of  the  exact  ABCs  is  often  translated  into  cumber¬ 
someness  and  high  computational  cost.  Thus,  along  with  the  accompanying  geometric  restrictions,  it  may  be 
regarded  as  a  serious  limitation.  The  alternative  is  provided  by  various  approximate  local  methods,  which 
typically  meet  the  other  usual  requirements  of  the  ABCs  besides  minimization  of  the  error  associated  with 
the  domain  truncation.  These  other  requirements  are  low  computational  cost,  geometric  universality,  and 
implementation  without  difficulties.  Still,  the  basic  trend  in  terms  of  accuracy  remains  the  following:  Higher 
accuracy  for  the  boundary  procedure  requires  more  of  the  nonlocal  nature  of  exact  ABCs  to  be  somehow 
taken  into  account. 

In  fact,  almost  any  numerical  algorithm  for  setting  the  ABCs  can  be  thought  of  as  a  compromise  between 
the  two  foregoing  groups  of  requirements  that  in  a  certain  sense  contradict  one  another.  Shifting  the  balance 
towards  locality  and  practical  efficacy  often  implies  insufficient  accuracy;  shifting  it  to  the  other  end,  towards 
highly  accurate  nonlocal  techniques,  may  often  yield  cumbersome  and  all  but  impractical  algorithms.  It  is 
not  surprising,  therefore,  that  the  treatment  of  external  boundaries  in  modern  production  computations 
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typically  follows  the  first,  local,  path.  In  computational  fluid  dynamics  (CFD),  for  example,  only  a  few 
ABCs’  methodologies  out  of  the  wide  variety  proposed  to  date  can  be  regarded  as  commonly  used  tools. 
All  of  them  are  either  based  on  essential  model  simplifications,  e.g.,  local  quasi-one-dimensional  treatment 
in  the  vicinity  of  the  artificial  boundary,  or  obtained  as  a  localization  of  some  nonlocal  ABCs.  To  meet  the 
overall  accuracy  requirements  when  using  such  simple  boundary  procedures,  one  often  has  to  use  excessively 
large  computational  domains. 

An  indepth  review  of  different  ABCs’  methodologies  that  have  been  published  in  the  literature  over 
the  recent  years  is  available  in  the  paper  by  Tsynkov  [21].  Besides  the  general  review,  in  this  paper  we 
focus  on  the  group  of  approaches  associated  with  the  generalized  potentials  of  Calderon’s  type  and  the 
difference  potentials  method  (DPM)  by  Ryaben’kii  [22-24].  The  application  of  the  DPM  provides  a  new 
and  very  powerful  vehicle  for  developing  ABCs  in  different  settings.  In  the  framework  of  the  DPM,  the 
boundary  conditions  are  obtained  using  the  equivalent  boundary  parameterization  of  the  entire  variety  of 
exterior  solutions;  the  latter  parameterization  is  built  with  the  help  of  the  so-called  generalized  Calderon’s 
boundary  projections.  The  DPM-based  boundary  conditions  are  usually  global.  However,  when  applied, 
e.g.,  to  solving  the  steady-state  external  problems  in  CFD,  they  combine  the  advantages  relevant  to  both 
global  and  local  methods.  In  other  words,  the  principal  gain  from  using  the  DPM  is  that  the  method  allows 
one  to  simultaneously  meet  the  high  accuracy  standards  of  the  ABCs  and  the  requirements  of  geometric 
universality  and  easiness  in  implementation.  In  addition  to  the  review  [21],  there  are  original  publications  on 
the  DPM-based  ABCs,  we  mention  here  only  those  closely  related  to  the  steady-state  compressible  viscous 
aerodynamics,  see  [25-33]. 

The  DPM-based  ABCs  have  been  implemented  along  with  the  NASA-developed  multigrid  flow  solvers  in 
both  two  and  three  space  dimensions.  The  investigated  flow  regimes  range  from  the  very  low  (incompressible 
limit)  to  transonic  speeds,  include  different  geometries,  laminar  and  turbulent  cases,  and  sometimes,  relatively 
complex  flow  phenomena,  like  shock-induced  separation  and  jets.  Compared  to  the  standard  local  boundary 
conditions,  the  DPM-based  ABCs  provide  for  much  better  accuracy  and  substantially  smaller  computational 
domains ,  significantly  faster  multigrid  convergence  [32],  and  improved  robustness  of  the  overall  numerical 
procedure.  We  particularly  emphasize  the  improved  combined  performance  of  global  ABCs  and  multigrid 
flow  solvers  as  similar  behavior  has  also  been  observed  by  other  authors  [34-38]. 

In  fact,  the  effect  of  convergence  acceleration  when  the  multigrid  algorithm  is  supplemented  by  global 
ABCs  [32]  was  a  part  of  the  reason  for  conducting  the  current  study.  Of  course,  there  are  significant 
differences.  In  our  previous  work  we  have  used  production  flow  solvers  that  by  themselves  have  room  for 
improving  the  performance,  at  least  from  the  theoretical  standpoint.  In  this  work,  we  use  the  solver  that  has 
already  been  optimized  for  performance  and  thus  its  convergence  rate  cannot  be  improved.  Therefore ,  we 
primarily  aim  at  demonstrating  the  possibility  to  significantly  reduce  the  size  of  the  computational  domain 
(and  consequently ,  the  grid  dimension)  while  maintaining  the  solution  accuracy,  as  well  as  optimal  multigrid 
convergence  rate. 

1.3.  Motivation  for  the  current  study.  As  outlined  above,  the  new  factorizable  schemes  and  cor¬ 
responding  multigrid  solvers,  as  well  as  global  DPM-based  ABCs,  have  independently  demonstrated  per¬ 
formance  superior  to  that  of  the  standard  methodologies  in  various  settings.  Besides,  in  many  cases  the 
application  of  global  ABCs  has  helped  to  speed  up  the  multigrid  convergence.  Therefore,  it  seems  most  nat¬ 
ural  to  try  and  combine  the  two  techniques  —  factorizable  discretizations  with  advanced  multigrid  and  global 
ABCs  —  with  the  hope  of  obtaining  a  joint  methodology  that  would  on  one  hand  be  capable  of  producing 
accurate  flow  solutions  on  domains  of  substantially  reduced  size  (compared  to  what  the  standard  methods 
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allow,  see,  e.g.,  [26,31])  and  on  the  other  hand,  will  not  hamper  the  unimprovable  multigrid  convergence 
rate  displayed  by  the  new  type  of  solvers  [14-16]. 

1.4.  Specific  objective.  Using  the  simplest  existing  versions  of  both  the  factorizable  scheme  and 
nonlocal  boundary  conditions  as  a  testing  ground,  we  intend  to  thoroughly  work  out  all  the  details  of 
merging  the  two  methodologies  (from  theoretical  issues  to  the  implementation)  and  then  experimentally 
study  the  overall  performance  for  a  series  of  simple  problems  that  allow  for  a  direct  comparison  with  the 
exact  solution.  We  emphasize  that  “the  simplest  existing  versions”  does  not  always  mean  the  latest  most 
universal  and  advanced  ones,  but  rather  those  that  have  been  tested  and  that  most  easily  lend  themselves  to 
the  analysis  by  analytical  means.  Specifically,  we  employ  the  so-called  pressure-Poisson  formulation  of  the 
factorizable  scheme  [14-16]  for  calculating  the  two-dimensional  inviscid  incompressible  fluid  flow  around  an 
airfoil.  Analytical  solutions  for  such  flows  are  readily  available  through  the  use  of  the  conformal  mapping 
technique.  The  global  exact  artificial  boundary  conditions  are  also  constructed  semi-analytically  by  means 
of  the  conformal  mapping  and  Fourier  transform  along  the  boundary,  with  no  explicit  use  of  the  DPM.  The 
DPM  itself  will  be  used  at  a  later  stage,  when  the  boundary  conditions  of  the  same  quality  are  required  for 
more  complex  settings. 

2.  Description  of  the  algorithm. 


2.1.  Discretization  of  the  equations  and  multigrid  algorithm.  The  algorithm  used  here  is  based 
on  a  special  formulation  of  the  flow  equations  that  involves  the  so-called  Pressure  Poisson  Equation  (PPE). 
The  initial  version  of  the  algorithm  was  presented  in  [13].  The  PPE  was  derived  on  the  differential  level  and 
then  discretized.  So  were  the  corresponding  boundary  conditions.  It  was  realized  later  that  deriving  the 
discrete  PPE  directly  can  simplify  greatly  the  treatment  of  the  boundary  conditions  and  also  introduce  some 
desirable  features  into  the  scheme  (like  a  certain  conservation  property  that  appeared  to  be  crucial  for  the 
purpose  of  the  current  work,  see  Section  2.2.1).  This  algorithm  is  described  in  detail  for  both  unstructured 
triangular  and  structured  quadrilateral  grids  in  [15].  Computational  results  obtained  with  this  algorithm 
have  been  presented  previously  in  [14,16].  In  the  current  work,  only  structured  grids  have  been  used.  Since 
some  features  of  the  discretization  appeared  to  be  critical  for  the  purpose  of  this  work,  namely,  constructing 
the  global  ABCs  that  allow  to  maintain  the  optimal  convergence  rates  of  the  algorithm,  we  shall  describe  the 
discretization  scheme  and  the  algorithm  in  detail.  The  discretization  is  based  on  the  incompressible  Euler 
equations  in  primitive  variables: 


du  du  dp 

Udi  +  Vdi  +  Tx-"' 

dv  dv  dp 

ute+vdi  +  di~0' 


du  du 
dx  +  dy  ’ 


where  u  and  v  are  Cartesian  components  of  the  velocity  vector  u  in  the  x  and  y  directions,  respectively, 
and  p  is  the  pressure.  The  density  is  taken  to  be  one.  The  Euler  equations  can  be  written  in  the  vector  form 
as  follows: 


( u  *  grad)«  +  gradp  =  0, 
div  u  =  0. 
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(2.1a) 

(2.1b) 


We  shall  now  derive  the  PPE  formulation  of  the  equations.  We  subtract  from  the  momentum  equation 
(2.1a)  the  velocity  vector  times  the  continuity  equation  (2.1b)  and  apply  the  divergence  operator  to  the 
result: 


div  (gradp  +  ( u  •  grad)u  —  u  div  it)  ==  0. 

The  previous  equation  obviously  reduces  to 

A p  =  f  =  -div((u  -  grad) it  -  u  div  «).  (2.2) 

The  PPE  formulation  of  the  Euler  equations  is  given  by  replacing  the  continuity  equation  (2.1b)  by  equation 
(2.2),  i.e.,  by  using  the  system  (2.1a),  (2.2)  instead  of  system  (2.1a),  (2.1b);  this  replacement  is  obviously 
equivalent. 

The  reason  for  introducing  the  new  PPE  formulation  is  that  once  we  formally  consider  the  transport 
coefficients  in  the  momentum  equation  constant,  i.e.,  use  the  equation 

(a  ■  grad)u  +  gradp  =  0,  (2.3) 


where  a  =  const,  instead  of  (2.1a),  then  equation  (2.2)  becomes  a  Laplace  equation 

A p  =  —div  ((a  •  grad)u  —  a  div  u)  =  0 


(2.4) 


and  decouples  from  the  rest  of  the  system.  The  momentum  equation  (2.3)  can  be  looked  at  as  a  standard 
advection  equation  with  known  forcing  function  (pressure  gradient).  We  emphasize  that  in  this  constant- 
coefficient  case  the  right-hand  side  of  equation  (2.4),  i.e.,  the  term  —div  ((a  •  grad)it  —  a  div  w),  is  equal  to 
zero  no  matter  what  the  actual  value  of  div  u  is.  This  property  is  crucial  for  solving  discretizations  because 
on  the  discrete  level  the  value  of  div  u  is  typically  of  the  order  of  truncation  error  and  not  exactly  equal 
to  zero;  moreover,  when  solving  a  discretization  by  an  iterative  scheme,  initially  the  value  of  div  u  may 
simply  be  far  away  from  zero.  Note  that  in  the  general  nonlinear  case,  for  slowly  varying  velocity  fields  u 
(as  opposed  to  constant  a),  the  terms  on  the  right-hand  side  of  equation  (2.2)  can  still  be  considered  small 
(more  precisely,  of  a  lower  order  compared  to  the  Laplacian)  and  thus  regarded  as  subprincipal.  Therefore, 
the  coupling  they  introduce  is  weak,  and  consequently,  equation  (2.2)  can  be  considered  decoupled  from  the 
rest  of  the  system  for  the  purpose  of  constructing  the  relaxation  procedure. 

We  shall  now  discuss  the  issue  of  the  boundary  conditions  for  the  pressure.  The  order  of  the  PPE 
system  (2.1a),  (2.2)  is  higher  (by  one)  than  that  of  the  original  system  (2.1a),  (2.1b).  Therefore,  to  ensure 
the  well-posedness  of  the  boundary- value  problem,  an  extra  boundary  condition  for  the  pressure  is  required. 
However,  the  new  problem  should  still  be  equivalent  to  the  original  one.  This  implies  that  this  additional 
boundary  condition  needs  to  be  derived  from  the  boundary  conditions  specified  for  the  original  problem  and, 
possibly,  differential  equations  of  the  original  system. 

Consider  a  domain  H  with  the  boundary  Integrating  equation  (2.2)  over  Q  and  applying  Gauss’ 
theorem,  we  obtain 


JJ  div  (gradp  -f  ( u  •  grad ) u  —  u  div  u )  dx  dy  — 
n 

j>  (gradp  +  (u  •  grad)u  —  udiv  u)nds  =  0, 


(2.5) 
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where  n  is  a  unit  vector  normal  to  the  boundary  and  s  is  the  arc  length  of  dCl.  A  sufficient  condition  for 
equality  (2.5)  to  hold  is  that  at  the  every  point  (x,y)  6  dtt  the  following  is  true: 

(gradp  +  ( u  ■  grad ) u  —  u  div  n)  =  0.  (2.6) 


We  emphasize  that  as  we  essentially  require  in  (2.6)  that  the  overall  normal  flux  be  zero  at  every  point  on 
90,  then  we  do  not  need  to  distinguish  ahead  of  time  between  the  outward  and  inward  normals.  In  fact,  it 
will  be  convenient  to  assume  that  the  normal  is  always  directed  toward  the  center  of  curvature  of  90,  this 
will  allow  us  to  use  the  Frenet  formulae  (see  below) . 

Introducing  a  local  orthonormal  coordinate  frame  (n,s),  where  s  is  a  unit  vector  tangent  to  90,  and 
the  corresponding  components  un  and  ua  of  the  velocity  vector,  we  can  rewrite  equation  (2.6)  as  follows 


dp 

dn 


((«  •  grad  )ti)  +  undiv  u 


—  {u  •  (u  •  grad)n)  -I-  un div  u. 


(2.7) 


Relation  (2.7)  specifies  a  general  Neumann-type  boundary  condition  for  the  pressure.  If,  for  example,  a 
particular  portion  of  the  boundary  under  consideration  90'  is  a  solid  wail,  then  from  the  so-called  tangency 
boundary  condition  for  velocity:  un\QQ,  =  0,  we  obtain 


^  =  (w  ■  (w -grad)n)  =  u,({u  •  grad)n)s 


n  ’ 


where  1Z  =  7£(s)  is  the  curvature  radius  of  90'.  Note,  to  obtain  the  last  equality  in  the  previous  chain 
we  have  used  the  Frenet  formulae  for  plane  curves.  Clearly,  the  physical  interpretation  of  the  boundary 
condition 

|£  =  -i#  (2-8) 

dn  11  v  ' 

is  that  the  normal  acceleration  of  the  fluid  particles  moving  tangentially  to  the  wall  is  directed  along  the 
normal  n,  i.e.,  toward  the  center  of  curvature,  and  equal  to  the  local,  i.e.,  instantaneous,  centripetal  ac¬ 
celeration.  This  is  because  the  momentum  equation  projected  onto  the  normal  direction  would  read  that 
the  acceleration  of  the  fluid  particles  is  equal  to  the  minus  normal  pressure  gradient.  For  a  straight  wall, 
boundary  condition  (2.8)  obviously  reduces  to  the  zero  Neumann  boundary  condition  for  the  pressure: 


Let  us  now  turn  to  building  the  discrete  scheme.  The  momentum  equations  are  discretized  using  a 
standard  first-order  upwind- difference  approximation  to  the  advection  operator  and  a  second-order  central- 
difference  approximation  to  the  pressure  gradient.  To  prevent  the  resulting  discretization  from  degeneration 
near  stagnation  points  or  across  streamlines  aligned  with  the  grid,  a  regularizing  artificial  dissipation  term 
is  used.  Considering  the  term  in  the  ^-momentum  equation,  the  difference  operator  is 

u____  =— (wm_|_i/2jj  —  Vm+l/2}j)(um+lyj  ~  umj) 

+  “b  iym~l/2}j){um,j  ~  1^*)? 
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where  the  superscript  h  denotes  the  discrete  approximation  to  the  corresponding  differential  operator,  the 
quantities  um±  1/2,j  are  average  velocity  components  given  by 

um+l/2,j  = 

1,  , 

1/2, j  2  \U™tj  l,j)» 

and  vm±ij2j  are  artificial  viscosity  coefficients.  Analogous  expressions  can  be  written  for  other  operators: 
and  The  artificial  viscosity  coefficients  vm±i/2j  are  defined  as 

l/m±l/2J  =  max  (/zAttm±i/2,j ,  |^m±l/2,j|)  j 


where 


Aum+ 1/2,  j  —  (wm+l,i  wm,j)> 

Altm_i/2)j  —  (^m,j  l,j)» 

and  \x  >  1/2  is  an  adjustable  coefficient,  taken  to  be  0.7  in  the  current  work.  This  form  of  the  artificial 
dissipation  was  presented  in  [16]. 


(m-lj+1) 

• - 


(mj+l) 


(m-l/2,j+l/2)  ;  (m+l/2j+l/2) 

north 


(m+l,j+l) 

V 


-0- 


west 


I 

i  •' 

L  m>j : 
soulth 


eait 


(m4l,j) 


(m-l/2,j-l/2)  j  (m+l/2,j-l/2) 

+  - —  > —  ~  + 
(m-lj-l)  (m,j-l)  (m+IJ-1) 


Fig.  2.1.  Comjm£a£icma/  grid  segment  and  a  control  volume. 

The  discretization  of  the  Pressure  Poisson  Equation  is  based  directly  on  the  integral  formula  (2.5),  where 
the  domain  fl  is  taken  to  be  the  control  volume  Amj  with  the  boundary  dAmj  (see  Figure  2.1): 

j>  (gradp+  (u  ■grad),u  -  udivu)nds  —  0.  (2.9) 

fj 

Breaking  the  boundary  of  the  control  volume  dAm>j  into  four  parts:  “south,”  “north,”  “east,”  and  “west” 
faces,  we  can  rewrite  (2.9)  in  the  following  form: 

^[(■^771+1/2, j  ~  +  (Gm,j+ 1/2  -  Gm,j- 1/2)]  =  0, 
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where  h  is  the  mesh  size,  ^±1/2,7  are  the  east  and  west  fluxes,  and  Gmj±  1/2  are  the  north  and  south  fluxes 
through  the  faces  of  the  control  volume  Amj.  The  fluxes,  in  their  own  turn,  are  evaluated  using  the  midpoint 
quadrature: 

=h[(Pm,j  ~  Pm—l,j )  +  ^m-l/2,j  (wm-l/2tj+l/2  ~~  um-l/2,j-l/2) 
-Wm-l/2,j(Vm-l/21j+l/2  “  ^m-l/2,j-l/2)] 

Gm,j-l/2  —h[(Pm,j  ~  PmJ-l)  +Wm)j_i/2(wm+i/2,j-i/2  “  «m-l/2,j-l/2) 

1/2  (^m-f  l/2,j— 1/2  1/2, j— 1/2  )]  j 

where  the  quantities  with  fractional  indices  are  obtained  using  either  linear  or  bilinear  interpolation: 


Um-l/2,j  *“  (^m,j  4*  Wyn-l,j)/2, 
um—l/2J—  \  ~  {um,j  +  Um-l,j  H"  Um,j-1  4*  Wm— l,j— 1)/4- 


The  fluxes  Fm+i/2}j  and  Gm,j+i/2  are  evaluated  similarly. 

We  now  describe  the  derivation  of  the  discrete  PPE  at  a  boundary  point  (m,0)  (see  Fig.2.2),  assuming 
that  no  boundary  condition  on  the  pressure  is  prescribed  there  ahead  of  time.  As  in  the  continuous  case, 
the  discrete  boundary  condition  for  the  pressure  will  follow  from  the  existing  physical  boundary  conditions 
and  the  discrete  equations  of  the  finite-difference  scheme  that  we  build. 


(m-1,1) 


(m,l) 


(m+1,1) 


(m-1/2,1/2)  !  (m+1/2,1/2) 

nohh 


(m-1,0) 


(m,0) 


(m+1,0) 


Fig.  2.2.  Boundary  control  volume. 


We  rewrite  equality  (2.5)  considering  the  control  volume  J3m)o  as  the  domain  shown  in  Figure  2.2: 

j)  (gradp  +  {u  •  grad)u  -  udiv  u)nds  =  0.  (2.10) 

dBmr  0 

Again,  as  in  the  case  of  an  internal  control  volume,  we  break  the  integral  (2.10)  into  four  parts  —  fluxes 
through  the  four  faces  of  the  control  volume  £m>0,  and  approximate  each  one  of  them  using  the  numerical 
quadrature 

2  (^m+l/2,0  “  Fm- 1/2, o)  +  K^m,l/2  ~  Gm, o)  =  0.  (2.11) 

However,  as  the  “south”  face  represents  now  a  solid  wall,  we  set  the  corresponding  flux  to  zero: 

Gm,0  =  0.  (2.12) 
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This  implies  (see  equation  (2.6))  imposing  the  appropriate  Neumann  boundary  condition  on  the  pressure 


h~dn  ”  “  Vm~ l/2,o)  +  ^m,oK+l/2,0  “  um- 1/2, o)- 

Expression  (2.13)  can  obviously  be  interpreted  as  a  finite-difference  approximation  of  the  equation 

dp  dv  du 

s;+“k-,’8;  =  0' 


(2.13) 


which  is  a  specific  version  of  equation  (2.6)  for  the  straight  solid  wall  boundary  y  =  0. 

The  flux  Gm)  1/2  in  (2.11)  is  defined  in  the  same  way  as  for  the  internal  control  volume  case.  For  the 
west  flux  we  use  the  following  formula 

Fm- 1/2,0  —  2  (Prny0  ”  Pm~  1,0 )  “  wm- 1/2,0  l/2,l/2  “  ^m-l/2,o) 

+^m-l/2,o(^m-l/2,l/2  “  ^m-l/2,o)j 


and  the  east  flux  Fm+ i/2,o  is  evaluated  similarly  to  the  above. 

This  discretization  is  extended  to  the  body-fitted  grids  (generalized  coordinates)  in  the  same  standard 
way  as  any  control-volume  discretization.  Therefore,  we  do  not  elaborate  on  this  issue  here. 

2.2.  Construction  of  the  ABCs.  As  follows  from  the  previous  section,  there  are  three  flow  quantities 
to  be  determined  throughout  the  computation  —  pressure  p  and  two  velocity  components  u  and  v.  External 
boundary  condition  (i.e.,  closing  mechanism  for  the  discretized  equations)  will  be  needed  for  each  of  these 
three  quantities.  As  the  exact  solution  for  the  problem  under  study  (incompressible  inviscid  airfoil  flow)  is 
known,  we,  for  methodological  purposes,  will  consider  different  computational  strategies.  First,  we  will  be 
solving  for  the  pressure  only  while  keeping  the  velocity  field  on  the  grid  as  prescribed  by  the  exact  solution. 
Next,  we  will  be  solving  for  all  three  flow  quantities.  In  both  cases,  we  will  compare  the  results  obtained  with 
the  global  ABCs  against  those  obtained  with  the  exact  boundary  conditions  of  the  Dirichlet  type  (Dirichlet 
data  at  the  external  boundary  taken  from  the  exact  solution)  and  on  the  other  hand,  against  the  results 
obtained  with  the  standard  local  ABCs.  Accordingly,  we  construct  the  global  ABCs  separately  for  the 
pressure  and  for  the  velocities.  The  boundary  conditions  are  obtained  using  the  separation  of  variables  along 
the  boundary;  this  approach  is  well-known  and  can  be  regarded  standard  for  deriving  global  ABCs,  see  [21]. 

2.2.1.  Pressure  ABCs.  In  the  pressure-Poisson  formulation  that  we  are  using,  the  elliptic  factor  in 
the  system  is  the  Poisson  equation  for  the  pressure 


A p=f,  (2.14) 

which  is  to  be  solved  on  the  infinite  domain  exterior  to  the  airfoil  subject  to  the  condition  of  boundedness 
of  the  solution  at  infinity  and  the  Neumann  boundary  condition 


dp 

dn 


r 


(2.15) 


on  the  airfoil  surface  T.  As  mentioned  in  Section  2.1,  boundary  condition  (2.15)  is  simply  a  reformulation 
of  the  momentum  equation  in  the  direction  normal  to  the  solid  wall,  the  function  ip  represents  centripetal 
acceleration  (depends  on  the  tangential  velocity),  the  continuous  explicit  expression  for  %p  is  given  in  (2.8) ,  and 


n 


the  discretization  is  built  as  shown  in  (2.13).  The  only  difference  that  to  be  emphasized  between  boundary 
condition  (2.15)  and  considerations  of  Section  2.1,  see  formula  (2.8),  is  that  formerly  we  have  considered 
the  normal  n  always  pointing  to  the  center  of  curvature  of  T,  whereas  n  in  equation  (2.15)  is  the  normal 
to  F  with  the  fixed  direction  toward  the  interior  of  the  airfoil  (i.e.,  external  normal  for  the  flow  domain). 
Therefore,  the  sign  of  the  centripetal  acceleration  in  the  function  ip  will  depend  on  whether  the  curve  T  is 
convex  or  concave  at  every  particular  location.  Similarly  to  i/>,  the  right-hand  side  /  of  equation  (2.14)  also 
depends  on  the  velocities,  the  explicit  expression  for  /  is  given  in  equation  (2.2).  Note,  a  more  general  and 
more  recent  formulation  of  the  scheme  involves  the  equation  for  the  velocity  potential  rather  than  pressure 
as  the  elliptic  factor.  For  incompressible  flows,  this  equation  is  always  homogeneous,  which,  in  particular, 
makes  the  construction  of  the  ABCs  conceptually  more  straightforward.  In  the  current  formulation,  as 
will  be  seen,  taking  care  of  the  inhomogeneity  /  requires  special  attention.  The  corresponding  experience, 
however,  is  going  to  be  useful  for  the  future  analysis  of  the  compressible  case  in  a  similar  framework. 

Let  us  now  denote  by  ft  the  unbounded  domain  on  which  we  are  solving  equation  (2.14),  then  F  =  dft. 
An  obvious  necessary  condition  for  solvability  of  problem  (2.14),  (2.15)  with  p  bounded  at  infinity  is  the 
equality  of  the  sum  of  all  sources  inside  the  domain  ft  to  the  total  flux  through  its  boundary  T,  i.e., 


q  r 


ipdT. 


(2.16) 


As  when  constructing  a  standard  analytic  solution  for  the  incompressible  airfoil  flow,  for  the  purpose  of 
constructing  the  ABCs  we  use  the  conformal  mapping  between  the  domain  S7  (exterior  of  the  airfoil)  of  the 
variable  2  =  x  4-  iy  6  C  and  the  exterior  of  the  unit  disk  {|£|  >  1 1,  (  =  £  +  ip  G  C}.  Here,  however,  we  will 
not  need  to  know  the  explicit  form  of  the  conformal  mapping,  it  will  be  sufficient  to  know  that  V  is  mapped 
onto  the  unit  circle  |£|  =  1  and  z  =  oo  corresponds  to  (  =  oo.  (For  uniqueness,  we  also  need  a  third  real 
parameter,  which  can,  for  example,  be  the  value  of  arg£'(oo),  which  corresponds  to  prescribing  the  angle 
of  rotation.)  It  is  also  important  that  the  Laplace  operator  of  (2.14)  does  not  change  with  the  conformal 
mapping.  As  concerns  the  notations,  we  will  for  simplicity  keep  the  same  symbols  p  and  /  for  the  functions  in 
the  new  coordinates  (£,77),  always  assuming  however  that  p  =  p(x(£,r)),y{£,r)))  and  /  =  /(#(£,  ??),  2/ (£,??))• 
Note,  the  actual  transform  that  we  use  is  given  by  formula  (3.1)  in  Section  3.1,  this  conformal  mapping 
maps  the  exterior  of  the  airfoil  onto  the  exterior  of  the  disk  of  radius  1.1  centered  at  (—0.1, 0)  (instead  of  the 
exterior  of  the  unit  circle).  This  is  only  a  technical  difference  that  will  not  affect  the  forthcoming  analysis 
and  conclusions. 

On  the  complex  plane  of  the  variable  £  =  £  -1-  ip  we  introduce  polar  coordinates  (/>,  6)  so  that  equation 
(2.14)  and  boundary  condition  (2.15)  transform  into 


1  9  /  dp\ 

pdp  \P dp) 


(2.17) 


and 

=  m,  (2-i8) 

p~  1 

respectively.  The  external  artificial  boundary  on  the  plane  £  will  be  a  circle  of  radius  R  centered  at  the 
origin:  p  =  R.  As  we  will  see  in  Section  3.1,  the  grid  used  for  the  numerical  integration  of  the  flow  equations 
is  obtained  by  transforming  a  polar  grid  from  the  plane  £  back  to  the  plane  2.  This  transform  maps  the 


dp 

dp 
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circle  p  =  R  onto  the  outer  coordinate  line  of  the  curvilinear  boundary-fitted  conformal  grid  built  around 
the  airfoil.  On  the  circle  p  =  R,  we  need  to  construct  the  ABCs  that  would  guarantee  that  the  solution 
of  equation  (2.17)  with  boundary  condition  (2.18)  found  for  1  <  p  <  R  can  be  smoothly  extended  beyond 
p  =  R  to  the  entire  infinite  domain  p  >  1  so  that  the  extended  solution  satisfy  the  same  differential  equation 
and  wall  boundary  condition  and  also  be  bounded  at  infinity. 

When  constructing  the  ABCs,  we  will  start  with  considering  a  particular  case  of  f(p,6 )  =  0  for  p  >  R. 
The  motivation  for  that  is  twofold.  First,  it  will  always  be  the  case  for  the  new  formulation  of  the  scheme 
for  incompressible  flows.  Second,  in  many  other  situations,  the  right-hand  side  /,  although  not  exactly 
zero,  decays  sufficiently  fast  in  the  far  field  so  that  by  neglecting  it  outside  the  computational  domain  (i.e., 
outside  the  artificial  boundary)  one  does  not  introduce  large  errors.  Experimentally,  it  is  always  possible  to 
see  whether  or  not  the  assumption  of  homogeneity  in  the  far  field  is  acceptable  [26,27,30,31,33]. 

Basically,  the  case  when  the  governing  equations  are  homogeneous  in  the  far  field  is  central  for  most 
ABCs’  methodologies,  see  [21].  For  the  particular  formulation  under  study,  having  built  the  homogeneous 
ABCs  we  will  also  show  how  one  can  explicitly  take  into  account  the  inhomogeneity  /,  see  (2.17). 

The  homogeneous  boundary  conditions  are  most  easily  obtained  in  the  continuous  formulation.  We 
Fourier  transform  equation  (2.17)  in  9  and  as  f(p,  9)  =  0  for  p  >  R,  arrive  at  the  following  family  of  ordinary 
differential  equations 


l_d 

pdp 


p>R, 


k  =  0,  ±1,  ±2, . . .  , 


(2.19) 


parameterized  by  the  wavenumber  k.  For  k  ^  0,  the  corresponding  homogeneous  equation  from  (2.19)  has 
two  linearly  independent  solutions:  p^  =  pw  and  p[2)  =  p~w.  Boundedness  at  infinity  implies  that  all 
modes  p^]  for  all  k  =  ±1,  ±2, . . .  should  be  excluded,  which  is  equivalent  to  the  following  countable  set  of 
conditions  written  in  terms  of  the  Wronskians: 


Pk 

_ 1 

dpk 

42) 

=  0, 

k  =  ±1,±2,...  . 

(2.20) 

dp 

dp 

Conditions  (2.20)  can,  in  fact,  be  imposed  at  any  location  p,  at  which  equation  (2.17)  is  homogeneous;  setting 
(2.20)  at  p  =  R  immediately  yields 


dpk 

dp 


I  p=R 


\p=R 


=  0, 


k  =  ±  1,±2,...  . 


(2.21a) 


As  concerns  k  —  0,  the  general  solution  of  the  corresponding  homogeneous  equation  from  (2.19)  is  C\  +C2  In  p- 
Boundedness  at  infinity  implies  C2  =  0,  which  leads  to  the  following  boundary  condition 


=  0 ,  (2.21b) 

dp  p=R 

which  can  also  be  formally  obtained  from  (2.21a)  by  letting  k  =  0.  The  countable  sequence  of  relations 
(2.21b),  (2.21a)  for  k  =  0,  ±1,  ±2, . . .  constitutes  the  full  set  of  exact  ABCs  at  the  external  artificial  boundary 
p  =  R.  Boundary  conditions  (2.21)  are  applicable  provided  that  f(p,0)  =  0  for  p>  R.  They  guarantee  that 
the  solution  calculated  for  p  <  R  can  be  smoothly  and  uniquely  complemented  to  the  entire  infinite  domain 
so  that  the  complement  solve  equation  (2.17)  and  be  bounded  at  infinity. 
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Note,  after  Fourier  transforming  back  to  the  physical  space,  the  ABCs  obtained  by  the  separation  of 
variables  along  the  boundary  typically  become  nonlocal  (e.g.,  ABCs  (2.21)).  Most  often,  these  ABCs  (in  the 
physical  space)  are  represented  by  means  of  pseudodifferential  operators,  see  [21].  For  the  particular  case  of 

(2.21),  however,  a  representation  via  singular  integrals  is  also  available;  it  was  obtained  by  Loncaric  in  [39] 

0 

and  involves  the  kernel  In  sin  -  . 

Let  us  now  consider  the  inhomogeneous  case,  when,  generally  speaking,  /(p,  6)  ^  0  for  p  >  R .  After  the 
Fourier  transform  in  9  we  obtain  (cf.  (2.19)) 

Id  (A)  -  = /*(,),  >!  fc  =  0, ±1, ±2, . . .  .  (2.22) 

p  dp  V  dp  j  p 2 

To  discuss  the  solvability  issues,  we  will  also  need  to  Fourier  transform  (2.18),  which  yields 

=$k,  fc  =  0,±l,±2,...  . 

dP  P= i 

We  start  with  the  analysis  of  the  case  k  =  0,  which  differs  from  the  analysis  for  all  other  k’s 
space  on  the  new  plane  £,  the  solvability  condition  (2.16)  obviously  transforms  into 

+°o 

/  fo{p)pdp  =  ^  =  i>o  ■  (2-24) 

We  now  integrate  equation  (2.22)  for  k  —  0  from  p  =  1  to  p  =  R  and  obtain: 


(2.23) 
.  In  Fourier 


dp  V  dp 


dp  =  p 


dp  \p=r  dp  \p=1 


=  /  fo(p)pdp  - 


Combining  formulae  (2.24)  and  (2.25)  immediately  yields 


dp  \o=i 


IX 

=  J  fo{p)pdp-i>o- 


Inhomogeneous  relation  (2.26)  will  replace  homogeneous  relation  (2.21b)  for  k  =  0  in  the  countable  sequence 
of  boundary  conditions  in  Fourier  space.  Note,  if  we  assume  for  a  moment,  as  before,  that  fo(p)  =  0  for 

/R  ^  r+oo 

fo{p)pdp=  I  fo(p)pdp,  and  because  of  (2.24)  the  right-hand  side  of  (2.26) 

vanishes  and  condition  (2.26)  transforms  back  into  (2.21b). 

Let  us  now  emphasize  a  very  important  circumstance.  Specifying  Neumann  boundary  conditions  for  k  = 
0  on  both  edges  of  the  interval  (i.e.,  for  p  =  1  and  p  =  R)  may,  generally  speaking,  be  problematic  from  the 
standpoint  of  solvability.  Indeed,  in  the  fully  homogeneous  case  neither  of  these  boundary  conditions  admits 
the  logarithmic  mode  and  both  admit  the  constant  mode.  Thus,  p0  =  const  will  solve  the  homogeneous 
problem 

1  d  f  dp0\  . 


=  0.  f- 

dp  p=R 
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Consequently,  the  corresponding  nonhomogeneous  problem  will  not  be  solvable  for  any  right-hand  side,  and 
in  case  it  is  solvable,  the  solution  will  not  be  unique.  The  type  of  the  boundary  conditions  (Neumann)  for 
k  =  0  cannot  be  changed/relaxed  because  at  both  p  =  1  and  p  =  R  these  boundary  conditions  are  obtained 
from  the  physical  considerations.  Therefore,  some  additional  solvability  conditions  are  introduced  that  limit 
the  general  admissible  scope  of  the  problem  data. 

As  we  have  seen,  for  the  original  problem  on  the  semi-infinite  line  p  >  1  the  solvability  condition  is 
required  specifically  because  the  boundary  condition  on  the  solid  wall  has  Neumann’s  type.  This  solvability 
condition  is  given  by  (2.24);  relation  (2.24)  actually  restricts  the  class  of  admissible  right-hand  sides  fo 
provided  that  ipo  is  given.  The  ABCs  set  at  a  finite  location  p  —  R  are  not  independent  from  the  interior 
problem;  it  is  rather  the  opposite  —  the  problem  with  these  ABCs  fully  inherits  the  properties  of  solvabil¬ 
ity/insolvability  of  the  original  problem.  As  has  been  shown  for  the  problem  with  the  ABCs,  the  solvability 
condition  (2.24)  is  rewritten  in  the  form  (2.26),  which  in  the  homogeneous  case  simplifies  to  (2.21b).  A 
convenient  circumstance  is  that  in  the  formulation  that  involves  the  ABCs  the  solvability  condition  (2.26) 
(or  (2.21b))  does  not  complement  all  other  conditions  of  the  problem,  it  is  rather  incorporated  naturally  (as 
a  particular  relation  in  the  Fourier  space  for  k  =  0)  into  the  countable  family  of  boundary  conditions  that 
we  build  in  Fourier  space  for  k  —  0,  ±1  =fc  2, . . . . 

Clearly,  the  homogeneous  boundary  condition  (2.21b)  can  be  interpreted  as  zero  flux  through  the  outer 
boundary  and  thus  it  corroborates  the  natural  physical  meaning  of  the  solvability  condition.  When  a  non-zero 
right-hand  /  side  extends  beyond  p  -  R,  we  replace  (2.21b)  with  the  nonhomogeneous  boundary  condition 
(2.26)  to  ensure  the  solvability.  The  meaning  of  (2.26)  is,  of  course,  the  same  —  it  is  conservation,  and  the 
difference  compared  to  (2.21b)  is  that  now  because  of  a  different  right-hand  side  the  flux  through  the  outer 
boundary  should  no  longer  be  zero.  The  important  thing  is  that  the  inhomogeneity  on  the  right-hand  side 
of  equality  (2.26): 


f  .  „  f+°°  „ 

J  fo(p)pdp  =  JR  fo(p)PdP 


(2.27) 


cannot  be  neglected  even  if  f0(p)  is  small  for  p  >  R,  see  (2.27).  As  has  been  mentioned,  the  latter  is  often 
the  case,  but  disregarding  the  aforementioned  inhomogeneity  will  make  this  problem  unsol vable  rather  than 
simply  introduce  a  small  error  into  the  boundary  conditions  and  thus  into  the  solution.  In  many  other 
situations  (see,  e.g.,  below)  introducing  such  errors  is  not  dangerous  from  the  standpoint  of  solvability,  and 
the  extent  of  the  corresponding  solution  deterioration  can  always  be  assessed  aposteriori  (see  [21]).  In  the 
particular  case  under  study,  strict  conservation  of  type  (2.26)  has  to  be  enforced  not  only  in  the  continuous 
framework,  but  also  on  the  grid;  this  issue  will  be  discussed  later. 

For  all  other  modes  except  k  =  0,  i.e.,  k  =  ±1,  ±2, . . . ,  it  is  easy  to  see  that  the  homogeneous  problem 


ld_  (  dpk\ 

p  dp  V  dp  ) 


k2 

~^Pk  =  0 

P 2 


dpk 

dp 


=  0, 


p= i 


dpk 

dp 


p=R 


0, 


p=R 


has  only  trivial  solution  (indeed,  both  p^  and  p^  are  eliminated  by  the  boundary  conditions)  and  therefore 
the  corresponding  nonhomogeneous  problem  is  solvable  for  any  right-hand  side  (inhomogeneities,  i.e.,  right- 
hand  sides  can  be  introduced  into  the  boundary  conditions  themselves  as  well).  Therefore,  the  situation 
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here  is  entirely  different  compared  to  the  case  k  =  0.  Namely,  boundary  conditions  (2.21a)  are  exact  when 
the  far  field  is  homogeneous,  f(p,0)  =  0  for  p  >  R.  If  a  far-field  inhomogeneity  is  introduced,  its  effect 
can,  in  principle,  be  incorporated  into  these  boundary  conditions  making  them  nonhomogeneous  as  well. 
However,  when  the  far-field  inhomogeneity  is  small  it  can  be  neglected,  which  will  not  give  rise  to  any 
solvability  concerns  for  k  ^  0,  and  the  corresponding  error  in  the  numerical  solution  will  then  be  estimated 
by  aposteriori  numerical  checks,  see,  e.g.,  [26,27,30,31,33]. 

Note,  the  incorporation  of  inhomogeneity  into  the  boundary  conditions  (2.21a),  k  ^  0,  is  not  going 
to  be  as  easy  as  obtaining  relation  (2.26)  instead  of  (2.21b).  A  general  “conceptual”  recipe  for  such  an 
incorporation  can  be  found,  e.g.,  in  [23];  but  it  basically  amounts  to  integrating  (essentially,  solving)  the 
nonhomogeneous  differential  equation  on  the  infinite  portion  of  the  domain  that  is  being  truncated  and 
replaced  by  the  ABCs.  This  is  exactly  what  we  are  trying  to  avoid;  besides,  it  is  not  feasible  at  any  rate 
unless  we  already  know  the  solution  or  at  least  the  closed  form  representation  for  f(p,6)  over  the  infinite 
domain. 

As  for  the  methodology  that  we  are  currently  describing,  we  neither  want  it  to  rely  heavily  on  the 
availability  of  the  exact  solution  nor  we  pursue  the  goal  of  assessing  the  performance  of  the  ABCs  in  the  regime 
when  their  exactness  is  “contaminated”  by  some  small  far-field  inhomogeneities  that  are  being  neglected.  As 
has  already  been  mentioned,  we  rather  aim  at  demonstrating  that  the  optimal  convergence  of  the  solver  can  be 
recovered  with  the  global  ABCs  in  the  simplest  (i.e.,  ideal)  regime.  Therefore,  in  the  numerical  experiments 
below  we  study  the  non-lifting  flows  so  that  the  entire  inhomogeneity  of  the  problem  is  concentrated  on  the 
zeroth  Fourier  mode  k  —  0  only.  Even  for  a  simple  lifting  airfoil  flow,  we  would  have  already  had  to  either 
neglect  the  far-field  inhomogeneities  for  k  =  ±1  and  thus  study  the  overall  performance  of  the  algorithm 
affected  by  this  approximation  as  well,  or  alternatively,  use  the  explicit  form  of  the  exact  solution,  substitute 
it  into  the  left-hand  side  of  (2.21a)  for  k  =  ±  1  and  this  way  create  the  correct  inhomogeneity  for  the  boundary 
conditions.  In  this  paper,  we  are  doing  neither  of  the  above.  As  concerns  studying  the  effects  of  dropping 
the  far-field  inhomogeneities,  an  investigation  of  this  type  for  the  more  general  compressible  formulation  of 
the  scheme  is  likely  to  become  a  subject  of  our  future  work. 

To  derive  the  ABCs  for  the  finite- difference  scheme  we  follow  a  procedure  similar  to  the  one  we  used 
for  obtaining  the  continuous  boundary  conditions.  The  grid  on  the  plane  C  is  a  polar  grid  uniform  in  6  and 
stretched  in  p\  when  mapped  back  onto  the  plane  z  it  yields  a  curvilinear  O-type  grid  fitted  to  the  airfoil 
surface.  The  grid  has  J  cells  in  the  radial  direction  with  the  nodes  pj ,  j  =  0, . . .  ,  J,  so  that  p0  —  1  and 
pj  =  Ry  and  M  cells  in  the  circumferential  direction  with  the  constant  grid  size  A 0  —  2? t/M  and  nodes 
Qm  =  mAd,  m  =  0, . . .  ,  M;  due  to  periodicity  the  angular  directions  6o  =  0  and  6m  =  coincide. 

A  second  order  accurate  discretization  of  the  flow  equations  that  we  use  is  outlined  in  Section  2.1.  It 
is  a  finite-volume  discretization  performed  directly  on  the  curvilinear  body-fitted  grid  in  the  physical  plane 
z.  It  is  easy  to  check  however,  that  if  the  exact  same  finite- volume  approach  was  applied  to  the  polar  grid 
on  the  model  plane  f,  then  for  the  pressure  Poisson  equation  (2.17)  it  would  result  in  the  following  natural 
central-difference  discretization: 

11/^  Pm,j-\- 1  Pm 

Ti  V;+1/2  A pi+1/2 

1  Pm+lijf  ~  2Pm,j  +  Pm—lJ  _  , 
pj  A02  ~ TmJ' 

where  pj+ 1/2  =  (pj+i  +  pj)/ 2,  pj-1/2  =  (pj  +  Pj- 1)/2,  &Pj+ 1/2  =  Pj+i  -  Pj>  &Pj- 1/2  =  Pj  ~  Pj- 1.  and 


,j  PmJ  Pm,j—  1 


+ 


(2.28) 
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A pj  =  (Apj -1/2  +  Apj+1/2)/ 2.  Finite-difference  equation  (2.28)  will  be  used  for  deriving  the  discrete 
counterpart  to  the  nonhomogeneous  boundary  condition  (2.26). 

Along  with  the  fully  discrete  equation  (2.28)  we  will  also  consider  a  semi-discrete  form  of  the  homogeneous 
equation  in  the  far  field: 


Id  (df^\  +  l^Prn+l  2j£m  +  Pm-1  =  p  m  =  0, . . .  ,  M  -  1.  (2.29) 

pdp  \  dp  )  p 2  A  vl 

Introducing  the  direct  discrete  Fourier  transform 

E  fc=-y  +  1---.y.  (2.30a) 

m— 0 

and  inverse  discrete  Fourier  transform 

M/2 

pm=  £  PkeikmA6,  m  —  ,M  —  1,  (2.30b) 

k=~M/  2+1 

on  the  grid,  we  reduce  (2.29)  to 


1  d_  (  dpk\ 
p  dp  V  dp  ) 


ak  - 

~  -4 Pk 
P2 


0, 


Oil. 


,2  _ 


fcA0 

2 


(2.31) 


For  any  particular  fc,  equation  (2.31)  looks  very  much  like  (2.19)  except  that  k 2  is  replaced  by  To 
obtain  the  boundary  conditions,  we  introduce  a  row  of  ghost  nodes  with  j  =  J  -f  1,  pj+i  =  p j  +  Apj+1/2  = 
JF£+ A/9/+1/2,  and  first  write  down  the  system  of  equalities  similar  to  (2.21a)  but  at  the  location  p  —  pj+1/2  — 
R  +  Apj+1/2/2  rather  than  p  =  R  (k  ^  0): 


dpk 

dp 


.  Kl  * 

H - PA; 

p=pj+iii  Pj+1/2 


p~pj + 1/2 


M  M 

=  0,  ife  =  -y  +  1,...  1,  1,...  ,  y.  (2.32) 


Then  we  discretize  (2.32)  consistently  with  the  central-difference  discretization  (2.28),  which  yields 


Pk,J+l~Pk,J  ,  Wk\  Pk,J+l  +Pk,J  _  n  11 

A  „  '  „  o  —  U)  o'-1) - J  - >  o  ’ 

^pj+1/2  .  Pj+1/2  2  2  1 


or 


Pk,. 


/M  |  ^J+l/2  ^ 

py+1/2  \ 

\  2  Apj+1/2  / 

1  ^  2 

^PJ+1/2/ 

,  M  ,  ,  M 

*  =  -y  +  !,•••  ,-l,  !,•••  > y • 


(2.33) 


(2.34a) 


The  remaining  mode  k  =  0  is  treated  separately  as  before;  in  the  homogeneous  case  instead  of  (2.21b)  we 
have: 

Po,j+i  ~  Po,J  _  q 
&Pj+ 1/2 
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or 


Pq,j+i  =Po,J-  (2.34b) 

Note,  as  in  the  case  of  the  continuous  boundary  conditions  (2.21b)  and  (2.21a),  relation  (2.34b)  can  also 
be  formally  obtained  from  (2.34a)  by  substituting  ak  =  0  for  k  =  0.  In  the  physical  discrete  space  we  can 
rewrite  the  system  of  equalities  (2.34a),  (2.34b)  as  one  matrix  relation 

pj-fi  =  F~Miag  {/?*}*>./  =  Tpj ,  (2.35) 


pj  and  pj+i  are  the  M-component  vectors  that  contain  the  grid  values  of  p  for  the  last  row  of  nodes  j  —  J 
and  the  row  of  ghost  nodes  j  =  J  +  1,  respectively,  and  F~ 1  and  F  are  M  x  M  matrices  that  denote  the 
inverse  (2.30b)  and  direct  (2.30a)  discrete  Fourier  transforms,  respectively. 

Note,  unlike  in  many  cases  from  our  previous  work  (see,  e.g.,  [25-27,30,31,33,40]),  boundary  condition 
(2.35)  is  not  built  so  that  to  guarantee  extension  of  the  interior  solution  beyond  the  artificial  boundary  as 
an  exact  solution  of  the  homogeneous  finite-difference  equation  (equation  (2.28)  with  fmj  =  0).  This  would 
require  obtaining  the  exact  discrete  counterparts  to  the  two  modes  pj^  and  p^  in  polar  coordinates  on 
the  stretched  radial  grid,  which  is  not  an  obvious  task.  Thus,  by  considering  the  semi-discrete  equation 
(2.29),  deriving  the  corresponding  continuous  boundary  conditions  (2.32),  and  then  approximating  them 
with  finite  differences  (2.33)  consistently  with  (2.28),  we  come  here  as  close  as  we  can  to  building  the  true 
exact  finite-difference  ABCs.  However,  as  we  always  maintain  the  second  order  of  accuracy  for  our  discrete 
approximations,  both  for  the  differential  equation  itself  and  for  the  boundary  conditions,  we  do  not  expect 
that  this  approach  may  be  problematic  form  the  standpoint  of  final  accuracy.  Our  numerical  experiments 
corroborate  that  the  ABCs  that  we  construct  this  way  are  capable  of  producing  very  high  accuracy,  quite 
comparable  with  the  one  obtained  when  the  ABCs  were  built  based  on  the  exact  solution,  see  Section  3. 

Moreover,  it  is  obvious  that  for  small  A0  and  small  k  (long  waves)  we  have  a\  «  ft2,  see  formula 
(2.31).  Therefore,  the  difference  between  the  truly  continuous  boundary  conditions  (2.21a)  and  semi-discrete 
boundary  conditions  (2.34a)  is  not  going  to  be  large  unless  all  modes  including  high  frequencies  on  the  grid 
(i.e.,  large  fc’s)  are  well  represented  in  the  solution.  Experimentally  we  have  observed  that  replacing  \ak\ 
by  simply  \k\  in  relations  (2.34a)  indeed  does  not  cause  any  noticeable  changes  in  the  numerical  solution 
when  calculating  the  non-lifting  airfoil  flows.  (The  exact  solution  for  such  flows  does  not  contain  any  angular 
modes  beyond  \k\  =  2.) 

Let  us  point  out  some  easy-to-see  properties  of  the  matrix  T  of  (2.35).  First,  this  matrix  is  circulant: 


to  tM- 1  tM- 2  •  •  ’  t\ 

tl  to  tM- 1  •  •  •  t2 

T  —  t2  ti  to  ...  *3  (2.37) 

tM-l  tM- 2  ^M- 3  t0 
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This  follows  immediately  from  the  periodicity  in  circumferential  direction.  Moreover,  the  matrix  T  is 
symmetric  because  tj  —  tM^j  for  j  =  1, . . .  ,  M  —  1;  these  equalities  simply  reflect  the  fact  that  the  response 
to  a  point  source  on  the  circular  boundary  will  be  the  same  in  both  clockwise  and  counterclockwise  directions. 
Finally,  all  eigenvalues  fa  of  the  matrix  T,  see  (2.36),  are  non-negative  and  moreover,  0  <  fa  <  1  for 
k  =  -M/2  +  1, ...  ,  M/2.  Indeed,  as  mentioned  in  Section  3.1  below,  the  grids  that  we  actually  use  for 
our  computations  are  stretched  in  the  radial  direction  so  that  to  maintain  the  cell  aspect  ratio  equal  to 
one  throughout  the  entire  domain.  This,  in  fact,  means  that  pj+1/2/ Apj+1/2  =  1/A0,  and  recalling  the 
definition  of  ,  see  (2.31),  we  conclude  that  the  first  expression  in  brackets  in  the  definition  of  fa,  see  (2.36), 
is  always  non-positive.  Thus,  fa  >  0.  The  second  inequality  for  the  magnitude  of  fa  follows  directly  from 
formula  (2.36). 

The  aforementioned  properties  may  be  beneficial  from  the  standpoint  of  multigrid  analysis,  as  well  as  a 
certain  type  of  implementation,  see  Section  2.3.  However,  it  yet  remains  to  be  seen  whether  similar  properties 
can  be  established  for  more  general  cases. 

When  a  far-filed  inhomogeneity  is  present,  we  need,  as  has  been  shown,  special  treatment  for  the  zeroth 
Fourier  mode  k  =  0.  First,  we  write  down  the  one-dimensional  finite-difference  equation 


1  1  /  POij+i  -  Pqj 

Pj  &Pj  \  17+1/2  &Pj+ 1/2 


Pj~  1/2 


Poj  ~  Pad- 


Apj- 1/2 


rh 


fo,jy 


(2.38) 


which  one  can  obtain  directly  from  (2.28),  and  notice  that  p0j  =  const  is  a  solution  to  the  homogeneous 
counterpart  of  (2.38).  This  is  the  only  case  an  exact  solution  to  the  discrete  homogeneous  equation  is  easy  to 
find.  The  second  linearly  independent  solution  to  (2.38)  has  to  grow  for  large  j9 s  as  on  any  finite  domain  it 
should  approximate  In  p\  this  second  mode  is  prohibited  by  boundary  condition  (2.34b)  for  the  homogeneous 
case  when  =  0  for  j  >  J.  Let  us  now  perform  summation  by  j  from  j  =  1  to  j  —  J  on  both  sides  of 
(2.38),  this  discrete  operation  is  an  analogue  of  the  integration  (2.25): 


t  (<■* 

j= 1 v 


X  .  .  P0.3+1  -  Po,j  _  P0,j  - 


P0,J4-1  -  P0,J 


7-1/2  / 


(2.39) 


Api/2 


j- 1 


From  (2.39)  we  conclude  that  we  need  to  replace  (2.34b)  by  the  following  nonhomogeneous  boundary  condi¬ 
tion: 


Po,j+i  ~  Po,J 


J 


E 


Pj&Pjfoj  +  Pl/2 


Pop  —  Po,o 
Api/2 


(2.40) 


which  is  a  discrete  counterpart  to  (2.26). 

Let  us  now  note  that  fm,j  has  a  rather  special  structure  —  it  can  be  represented  using  the  divergence 
theorem  as  fluxes  through  the  boundary  of  a  control  volume  Amj  shown  on  Figure  2.1,  see  the  discussion 
in  Section  2.1.  On  the  polar  grid,  the  area  of  this  control  volume  is 

^(^+1/2  ~  lj-1/2)  =  2k pjApj 
M  M  ‘ 


Now,  applying  the  divergence  theorem,  we  obtain 


f  . -a  n(Ff 

J™,3  — 


m+l/2,j  FL-l/2,j )  +  (tg»  J+l/2 


.  (2.41) 
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Notice,  unlike  in  Section  2.1,  the  fluxes  i^±1/2  j  and  j±  1/2  *n  eQuation  (2.41)  have  the  superscript  /; 

this  emphasizes  that  they  take  into  account  only  the  contribution  of  the  actual  right-hand  side  /  of  equation 
(2.2).  This  means  that  instead  of  evaluating  integral  (2.9)  using  the  control  volume  approach,  we  evaluate 
a  similar  integral  with  gradp  taken  out  and  thus  only  the  velocity’s  contribution  left. 

Calculating  now  the  zeroth  Fourier  component  of  /mj-,  see  equation  (2.30a),  and  using  equation  (2.41) 
we  obtain 

1  M— 1  1  /  M-l  M— 1  \ 

f°’i  =  wE  ^  =  TTcrAO'  (  Pj+1/2  EE  G™,}+ 1/2  “  ^'-!/2  EE  Gli,j-l/2  )  * 

771=0  ^  ^  \  771=0  m=0  / 

Substituting  the  previous  expression  into  (2.40)  yields 


PO  ,j+i  ^ 


M 


2^  ' 

771=0 


T77l,  J+l/2 


M 


(2.42) 


771=0 


Clearly,  the  discussion  in  Section  2.1  and  more  precisely,  the  requirement  (2.12)  of  zero  flux  through  the 
boundary  that  leads  to  the  boundary  condition  (2.13)  on  the  solid  wall,  implies  that 


Therefore,  we  obtain 


P1/2 

M 


M-l 


EE  Gm,  1/2  -  Pl/2 

771=0 


fto,i  —  £0,0 
Apl/2 


„  /V  -  iVfl  —  1 

P0,J+1  —  P0,J  _  1 
&PJ+1/2  M 

Relation  (2.40)  and  its  particular  form  (2.43)  represent  the  exact  conservation  on  the  grid  that 
for  the  solvability  of  the  discrete  problem.  The  inhomogeneity  on  the  zeroth  Fourier  mode,  see 
also  be  included  into  the  boundary  condition  in  the  matrix  form,  which  will  then  read 


(2.43) 

is  required 
(2.43),  can 


Pj+i  =  Tpj  + 


Apj+1/2  7 
2ttPj+i/2  ,Jj 


(2.44) 


where  f.tj  is  the  total  flux  due  to  the  forcing  term  through  the  outer  faces  of  the  control  volumes  centered 
at  the  J-th  gridline. 

From  the  standpoint  of  implementation,  boundary  conditions  (2.35)  or  (2.44)  are  matrix- vector  relations 
that  connect  the  values  of  the  solution  p  on  the  last  row  of  grid  nodes  j  —  J  and  the  row  of  ghost  nodes 
j  —  J- j- 1 .  These  boundary  conditions  are  equivalent  to  the  possibility  of  smoothly  and  uniquely  complement 
the  solution  from  inside  the  computational  domain  to  its  infinite  exterior  so  that  the  extension  solve  the  entire 
original  problem.  The  unique  choice  of  such  boundary  conditions  (within  the  accuracy  of  finite-difference 
approximation)  along  with  the  existence  of  the  conformal  mapping  between  the  planes  z  and  £  imply  that  the 
discrete  ABCs  (2.35)  or  (2.44)  obtained  originally  for  the  plane  (  can  be  used  on  the  plane  0  with  no  changes, 
simply  as  the  equalities  connecting  the  values  of  the  solution  on  two  rows  of  the  grid  through  matrix- vector 
multiplication.  This  property  is  fully  corroborated  experimentally,  see  Section  3,  even  so  the  finite-difference 
equation  for  the  pressure  on  the  actual  body-fitted  grid  may  differ  from  equation  (2.28)  in  the  sense  that 
the  stencil  may  contain  more  than  five  nodes  for  those  cells  that  noticeably  differ  in  shape  from  rectangles. 
(In  contrast  to  (2.28),  for  strongly  skewed  cells  diagonal  nodes  may  be  present  with  small  coefficients.) 

We  also  note  that  even  provided  that  the  problem  is  solvable,  ensuring  the  uniqueness  of  the  discrete 
solution  requires  special  attention.  Indeed,  as  we  have  seen  an  arbitrary  constant  can  be  added  to  the 
solution  on  the  zeroth  Fourier  mode.  The  issue  of  guaranteeing  the  uniqueness  throughout  the  computation 
is  touched  upon  in  the  implementation  Section  2.3. 
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2.2.2.  ABCs  for  velocities.  It  is  known  that  for  the  particular  case  under  study  the  quantities  u  and 
— v  satisfy  the  Cauchy-Riemann  equations: 


du  _  dv 
dx  dy  ’ 
du  _  dv 
dy  dx 


(2.45a) 


Let  us  first  note  that  equations  (2.45a)  imply  that  both  u  and  v  separately  satisfy  the  Laplace  equation. 
Thus,  we  essentially  could  have  applied  the  same  homogeneous  boundary  condition  that  we  have  constructed 
for  p,  see  (2.35),  to  velocities  as  well.  However,  we  will  rather  construct  the  ABCs  for  velocities  directly  on 
the  basis  of  the  Cauchy-Riemann  system.  This  will  simplify  the  analysis  for  the  zeroth  Fourier  mode,  which 
was  “problematic”  in  Section  2.2.1. 

As  before,  we  will  construct  the  ABCs  as  vector  relations  on  the  model  plane  (  and  then  use  them 
with  no  changes  on  the  physical  plane  z.  Also  as  before,  we  will  retain  the  same  notations  u  and  v  for  the 
functions  in  the  new  coordinates  (£,77),  always  keeping  in  mind,  however,  that  u  —  u(®  (£>*?)  >2/ (£>*?))  and 
v  =  u(z(£,7?),2/(£,?7)).  On  the  model  plane  (  =  £  +  ip  the  quantities  u  and  -v  also  satisfy  the  Cauchy- 
Riemann  equations: 


du  _  dv 
d£  dr] 
du  _  dv 
dp  d£ 


(2.45b) 


The  boundary  conditions  for  velocities  will  actually  be  constructed  on  the  basis  of  equations  (2.45b). 
Let  us  introduce  polar  components  of  the  velocity  vector  u 


up  — u  cos  9  +  v  sin  0  , 
ue  =  —  u  sin  6  +  v  cos  0  , 

and  rewrite  the  Cauchy-Riemann  equations  (2.45b)  in  polar  coordinates 

d(pup)  _  du& 
dp  "  dO  ’ 

d(pup)  __  dup 
dp  dO 

Denote  q  —  [up^up]1  and  represent  (2.47)  as  a  matrix  equation 

djpq]  _  .dq_ 

dp  39  ’ 

where 


For  the  Fourier  coefficients  qk  =  qk(p)  we  obtain  from  (2.48): 

1  =  ikAqk  ,  k  =  0,  ±1,  ±2, . . .  . 
dp 


(2.46) 


(2.47) 


(2.48) 


(2.49) 


(2.50) 


21 


We  now  diagonalize  the  matrix  A  of  (2.49)  and  obtain: 


AK  =  KA, 


(2.51) 


where 


A  = 


i  0 
0  -i 


K  = 


1  1 

—i  i 


Introducing  sk  =  K  1qk  qk-Kskj  we  get  instead  of  (2.50): 

d(ph) 


dp 


=  ikAsk  ,  k  =  0,  ±1,  ±2, . . .  . 


Equation  (2.53)  immediately  yields 


sk  ~ 


c(1k)p-k-1 

cfV-1 


(2.52) 


(2.53) 


(2.54) 


where  c[k ^  and  are  constants  (different  for  different  fc’s).  Going  back  from  representation  (2.54)  to  the 
old  variable  qk ,  we  have  (the  matrix  K  is  given  by  (2.52)): 


(2.55) 


cWp-k-i+c(k)pk-i 

upk 

_  -iC{k)p-k~l  +  iCik)pk-1  _ 

Qk  =  Ksk  = 

Finally,  introducing  the  inverse  to  the  transformation  (2.46) 

u  =  up  cos  6  —  uq  sin  9  , 
v  =  Up  sin  6  4-  uo  cos  0  , 

and  substituting  the  expressions 

u  =  J2^eike  , 

k 

up  =  '^TuPkeik8  , 
k 

into  (2.56)  we  obtain 


v  =  '£/vkeike  , 
k 

ue  =  J2uokeik6  , 

k 


k 

rei(k+l)e  _  ei(k-l)B 

'ei(k+l)B  +e*(fc- 1)0 

which  implies 


1.  i .  1.  i  . 

uk  —  2UP><-i  2Uf)k~1  2Upk+1  2U^k+1 

1  „  i  *  1  . 


—  2  ^pfc-i  2^fc+1 

The  last  step  is  to  substitute  (2.55)  into  (2.57): 

nk^C[k-1)p-k  +  Cik+1)pk  , 

ik+1)P' 


vk  =  -iCf ' ~k)p-k  +  iCik+1)-k 


(2.56) 


(2.57) 


(2.58) 
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Representation  (2.58)  is  valid  for  all  k’s  including  k  =  0.  Contrary  to  the  general  solution  of  the  Laplace 
equation  for  k  —  0:  C\  +  C2  In  p,  see  Section  2.2.1,  we  do  not  have  a  growing  logarithmic  mode  in  represen¬ 
tation  (2.58)  for  k  =  0.  This  constitutes  a  principal  difference  compared  to  the  previously  analyzed  case  and 
also  gives  the  primary  reason  why  we  construct  the  ABCs  for  velocities  directly  from  the  Cauchy-Riemann 
system  (2.45b)  rather  than  first  reduce  it  to  two  separate  Laplace’s  equations  for  u  and  v. 

The  physical  boundary  conditions  for  velocities  at  infinity  are  obvious:  the  vector  field  u  has  to  approach 
its  free-stream  value  t*o  as  \z\  — >  00.  As  C(°°)  =  the  value  of  u  at  infinity  on  the  model  plane  (  is  also 
and  thus  the  boundary  condition  at  infinity  on  the  plane  (  is  the  same  as  on  the  plane  z:  u  — >  «o  as 
| C|  — >  00.  This,  in  particular,  implies  boundedness  of  the  vector  field  u  at  infinity.  Then,  from  formula 
(2.58)  we  conclude  that  we  need  to  require  that  C =  0  for  k  >  0  and  c[k~ ^  =  0  for  k  <  0;  for  k  =  0 
both  and  are  allowed.  Therefore,  for  k  ^  0  in  the  continuous  framework  we  obtain  analogously 

to  (2.21a): 


duk 

dp 


p—R 


+  HUk 


=  0, 


p=R 


dvk 

dp 


p—R 


,  1*1, 
+  RVk 


=  0, 


p—R 


and  in  the  discrete  framework,  analogously  to  (2.34a): 


k  =  ±1,  ±2, . . .  , 
k  —  ±1,  ±2, . . .  , 


Uh,J+ 1 


fjM  +  pj+m  A 

\  2  ^Pj+l/2/ 


PJ+ 1/2  \ 
^Pj+1/2/  ’ 


Vk,J+l 


Hk\  +  Pj+ 1/2  \ 

V  2  &PJ+ 1/2/ 


Pj+1/2  \ 
^Pj+1/2/  ’ 


-1,  1,.. 


M 
2  ’ 


1,.. 


M 

y 


(2.59) 


(2.60) 


The  only  difference  between  (2.60)  and  (2.34a)  is  that  in  (2.60)  we  have  |fc|  instead  of  |afc|.  As  has  been 
shown,  oik  in  (2.34a)  comes  from  the  central-difference  second-order  discretization  of  the  Laplacian  (2.28) 
that  we  do  not  use  for  velocities.  It  has  also  been  mentioned  that  for  small  fc’s  (long  waves)  the  difference 
between  using  |A;|  and  |a*;|  even  in  (2.34a)  is  not  noticeable. 

For  k  =  0,  the  difference  between  the  boundary  conditions  for  velocities  and  the  corresponding  boundary 
conditions  for  the  pressure  is  significant.  As  we  do  not  have  to  worry  about  cancelling  the  growing  logarithmic 
modes  for  velocities,  and  as  we  know  the  free-stream  value  of  the  velocity  vector  no  =  [U,  V]4,  we  simply 
have 


wo,j+i  —  U  , 
£o,j+i  =  V  . 


(2.61) 


The  meaning  of  boundary  conditions  (2.61)  is  transparent  —  the  correct  value  of  the  solution  at  infinity, 
i.e.,  the  correct  constant,  is  picked  up  on  the  zeroth  Fourier  mode,  all  other  Fourier  modes  in  the  solution 
vanish  at  infinity. 

Formulae  (2.60)  and  (2.61)  constitute  the  global  ABCs  for  velocities.  It  is  interesting  to  note  that  as  the 
momentum  equations  are  actually  integrated  by  means  of  downstream  marching,  the  information  provided 
by  boundary  conditions  (2.60),  (2.61)  is  used  only  on  the  inflow  portion  of  the  artificial  boundary.  As, 
however,  the  computations  show,  the  values  of  u  on  the  outflow  portion  of  the  boundary  obtained  with  the 
global  boundary  conditions  very  well  agree  with  the  actual  values  obtained  by  the  downstream  integration 
through  the  domain,  see  Section  3.4. 
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2.3.  Combined  formulation  and  implementation.  The  values  of  the  solution  at  the  ghost  points, 
i.e.,  on  the  (J  4*  l)~th  gridline,  are  obtained  by  applying  the  ABCs  of  Sections  2.2.1  and  2.2.2  to  the  current 
approximation  on  the  outer  boundary,  i.e.,  on  the  J-th  gridline.  This  is  followed  by  a  relaxation  sweep.  The 
aforementioned  two  steps  form  the  combined  relaxation  sweep. 

We  apply  a  standard  multigrid  algorithm  known  as  the  Pull  Approximation  Scheme  (FAS).  The  only 
“non-standard”  aspect  that  requires  special  attention  is  the  application  of  the  ABCs  for  the  pressure  on  the 
coarser  grids,  because  these  ABCs  involve  the  conservation  issue  outlined  in  Section  2.2.1. 

The  pressure  equation  on  the  finest  grid  can  be  written  as  follows 

Lhph  =  fh , 

where  Lh  is  the  discrete  Laplacian  and  fh  is  the  forcing  term  as  defined  earlier.  The  fine  grid  equation 
residual  is  given  by 

rh  =  }h  —  Lhph. 

Then,  the  FAS  coarse  grid  equation  is 

Lhph  = fH, 

where 

/fl  =  W+/fr», 

and  Iff  and  Iff  are  the  restriction  operators.  The  operator  Iff  is  a  standard  injection  operator ,  and  the 
operator  iff  is  also  a  standard  Fall- Weighting  operator. 

Let  us  note  that  since  iff  is  a  “conservative”  operator,  the  coarse-grid  right-hand  side  fH  possesses 
the  same  telescopic  property  as  the  fine-grid  right-hand  side  fh  has.  This  means  that  it  is  composed  of 
flux  contributions  through  the  control  volume  boundary,  and  for  a  pair  of  neighboring  control  volumes  the 
contributions  from  fluxes  through  the  common  interface  cancel  one  another.  Therefore,  by  summing  up  the 
source  terms  fH  on  the  coarse  grid  we  again  arrive  at  the  total  flux  through  the  outer  faces  of  the  control 
volumes  centered  at  the  nodes  of  the  J-th  gridline,  i.e.,  the  quantity  2m=o  ^m,j+ 1/2 »  see  f°rinula  (2.43), 
(the  index  J  refers  here  to  the  coarse  grid)  in  the  same  way  as  we  have  done  that  on  the  fine  grid. 

As  the  pressure  in  an  incompressible  flow  problem  is  defined  only  up  to  an  additive  constant,  the 
uniqueness  of  the  solution  is  achieved  by  keeping  its  value  fixed  and  equal  to  an  initially  prescribed  quantity 
at  one  of  the  nodes  on  the  coarsest  grid  throughout  the  entire  computation. 

We  should  also  mention  that  another  way  of  implementing  the  non-local  ABCs,  which  has  not  been 
actually  tried  yet,  consists  of  using  relation  (2.44)  to  actually  eliminate  the  ghost  cells  and  thus  solve 
iteratively  only  for  the  variables  on  the  original  grid.  The  properties  of  the  ABCs’  matrix  T  obtained 
in  Section  2.2.1,  see  the  discussion  right  after  formula  (2.37),  may  then  help  to  establish  and  analyze  the 
properties  of  the  overall  system  matrix  from  the  standpoint  of  convergence  and  convergence  rate  of  the 
(multigrid)  iterations. 

3.  Numerical  experiments. 

3.1.  Test  problems.  One  test  problem  considered  here  is  the  incompressible,  irrotational  flow  around 
a  symmetric  airfoil.  A  circular  cylinder  is  transformed  into  an  airfoil  by  the  Karman-Trefftz  transformation, 
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where  £  is  the  complex  coordinate  in  the  circle  plane,  z  is  the  complex  coordinate  in  the  physical  plane, 
and  e  is  the  trailing-edge  angle.  The  cylinder  in  the  £-plane  is  centered  at  C  =  (-0-1, 0)  with  a  radius  of  1.1. 
The  trailing  edge  angle  e  is  10°.  This  generates  an  airfoil  with  a  thickness  ratio  of  approximately  15%  in 
the  2-plane. 

The  grids  are  generated  by  using  a  uniform  azimuthal  spacing  around  the  cylinder  in  the  £-plane.  The 
spacing  in  the  radial  direction  is  stretched  such  that  the  cell  aspect  ratio  is  equal  to  one  everywhere.  The 
dimensions  of  the  grids  that  we  have  used  are  given  in  the  next  section,  see  Table  3.1. 

Besides  the  airfoil  flow,  and  primarily  for  the  reasons  of  comparing  the  numerical  performance  for 
different  geometries,  we  have  also  calculated  an  irrotational  flow  around  a  circular  cylinder.  This  is  the  same 
cylinder  as  the  one  we  substituted  into  the  transformation  (3.1)  to  obtain  the  Karman-Trefftz  airfoil.  Thus, 
for  the  cylinder  case  the  planes  2  and  £  coincide  and  the  flow  around  the  cylinder  is  calculated  on  the  actual 
polar  grid. 

Only  non-lifting  solutions  are  considered  hereafter  for  both  the  airfoil  and  cylinder,  so  the  freestream 
flow  is  aligned  with  the  direction. 

3.2.  Computational  setting.  The  actual  dimensions  and  approximate  geometric  sizes  of  the  grids 
that  we  have  used  for  our  calculations  are  presented  in  Table  3.1. 


Table  3.1 

Grid  dimensions  and  sizes 


Grid  dimension:  angular  x  radial 

129  x  129 

129  x  65 

129  x  33 

129  x  17 

129x9 

Approximate  location  of  the 
outer  boundary  for  the  cylinder 
case  measured  in  diameters 

272.5 

11.5 

2.5 

1.1 

0.75 

Approximate  location  of  the 
outer  boundary  for  the  airfoil 
case  measured  in  chords 

135 

6 

1.25 

0.55 

0.37 

Approximate  locations  of  the  boundary  in  Table  3.1  are  given  with  respect  to  the  center  of  the  cylinder  (i.e., 
circle)  or  airfoil;  in  the  latter  case  the  location  closest  to  the  airfoil  surface  is  presented. 

On  each  grid  shown  in  Table  3.1  we  have  calculated  a  non-lifting  incompressible  inviscid  flow  with  a  given 
free-stream  speed.  The  discrete  flow  solution  is  calculated  by  applying  multigrid  iterations  to  the  elliptic 
part  of  the  factorized  system  (pressure)  and  using  a  downstream  marching  algorithm  for  the  advection  part 
(velocities),  see  Section  2.3.  We  employ  a  full  approximation  scheme  (FAS)  multigrid  W(2,l)  cycle,  with 
full  weighting  restriction  and  full  coarsening  [10].  The  number  of  nested  grids  for  each  computation  is 
determined  by  the  dimension  of  the  finest  grid  (see  Table  3.1),  i.e.,  the  number  of  possible  subdivisions  by 
two.  On  the  coarsest  grid  we  find  the  “exact”  solution  by  performing  sufficiently  many  relaxation  sweeps. 
For  a  pure  Laplace  or  Poisson  equation  on  the  grids  of  the  kind  presented  in  Table  3.1,  the  smoothing  rate 
for  a  lexicographic  Gauss-Seidel  relaxation  is  0.5  per  one  sweep  [10].  Therefore,  for  a  W(2,l)  cycle  the 
“predicted”  or  ideal  convergence  rate  is  0.5<2+1)  —  0.125  per  cycle,  i.e.,  the  residual  should  drop  by  a  factor 
of  1/8  in  a  single  multigrid  cycle. 

Actually,  we  have  calculated  two  types  of  flow  solutions.  First,  we  were  solving  for  the  pressure  only 
while  keeping  the  velocity  field  frozen  on  the  grid  with  the  values  taken  from  the  exact  solution.  Next,  we 
solved  the  full  Euler  system  as  well.  The  ideal  convergence  rate  of  about  0.125  per  cycle  was  indeed  observed 
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in  our  computations  when  we  were  solving  for  the  pressure  only,  see  Section  3.3.  As  concerns  the  true  Euler 
solution,  if  we  were  to  fully  separate  the  elliptic  and  hyperbolic  factors  in  the  Euler  equations,  we  should  have 
seen  the  exact  same  convergence  rate.  In  fact,  we  observed  a  somewhat  slower  convergence  when  solving  the 
full  Euler  equations,  see  Section  3.4.  The  reason  for  the  discrepancy  is  not  related  to  the  treatment  of  the 
external  artificial  boundary  as  the  slowdown  takes  place  for  all  three  types  of  the  ABCs  that  we  employ,  see 
below.  Convergence  deterioration  that  we  observe  is  rather  caused  by  the  following.  It  has  been  mentioned  in 
Section  2.1  that  the  subprincipal  terms  in  this  equation  can  be  disregarded  when  constructing  the  relaxation 
scheme.  However,  near  the  stagnation  point  this  may  be  not  legitimate.  A  possible  implication  of  this,  i.e., 
loss  of  efficiency,  is  discussed  in  the  recent  paper  [16]. 

For  each  variant  of  the  computation  (determined  by  the  geometry  and  the  grid,  see  Table  3.1),  we  have 
used  three  types  of  external  ABCs.  The  first  boundary  conditions  are  of  the  Dirichlet  type  for  all  three 
flow  quantities,  the  actual  data  to  specify  at  the  outer  boundary  are  taken  from  the  available  exact  solution. 
This  is  apparently  the  best  possible  treatment  for  the  artificial  boundary  and  henceforth  these  ABCs  will 
be  referred  to  as  exact  The  second  boundary  conditions  are  globed  ABCs  described  in  Section  2.2,  hereafter 
they  will  be  referred  to  as  global  Finally,  the  third  boundary  conditions  are  local ;  they  are  set  as  follows. 

In  the  course  of  the  relaxation  procedure,  the  residuals  of  the  pressure  equation  are  evaluated  near  the 

dp 

outer  boundary  on  “halves-cells”  using  the  condition  —  =  0,  where  n  is  the  normal  direction  to  the 

outer 

boundary.  Besides,  at  each  point  on  the  outer  boundary  we  specify  the  flow  angle,  which  is  approximately 
taken  to  be  equal  to  that  at  infinity  (i.e.,  to  the  angle  of  attack).  Finally,  to  calculate  the  magnitude  of 
the  velocity,  we  use  the  conservation  of  the  total  pressure  and  the  corresponding  value  of  the  static  pressure 
that  has  just  been  updated.  The  local  approach  is  obviously  the  simplest  of  the  three  in  the  sense  that  local 
ABCs  are  easier  to  construct  and  implement  than  the  global  ones  and  at  the  same  time  they  can  be  used 
when  the  exact  solution  is  not  available. 

The  purpose  of  conducting  the  foregoing  series  of  computations  on  different  grids  (see  Table  3.1)  with 
different  ABCs’  is  to  compare  the  multigrid  convergence  rates  as  they  depend  on  the  type  of  ABCs  and 
domain  size.  We  will  also  be  comparing  the  accuracy  of  the  numerical  solutions  as  it  depends  on  the  domain 
size  and  the  type  of  the  ABCs.  These  two  characteristics  — -  convergence  rate  (i.e.,  numerical  efficacy)  and 
accuracy  —  are  obviously  of  the  foremost  importance  when  designing  any  numerical  algorithm. 

Since  the  solutions  that  we  have  calculated  develop  neither  lift  nor  drag,  the  accuracy  is  assessed  not 
by  examining  integral  characteristics  of  the  flow,  like  lift  and  drag  coefficients  (as  typical  for  many  CFD 
studies),  but  rather  by  comparing  the  actual  error  profiles  for  different  flow  quantities. 

3.3.  Solving  for  the  pressure  only.  In  this  section  we  describe  and  analyze  the  numerical  results 
for  two  cases  —  non-lifting  flow  past  a  circular  cylinder  and  non-lifting  flow  past  an  airfoil.  As  has  been 
mentioned,  the  velocity  field  on  the  grid  for  this  series  of  computations  is  considered  known,  we  actually 
take  it  from  the  exact  solution  obtained  using  classical  complex  variable  technique. 

3.3.1.  Circular  cylinder.  In  Figure  3.1  we  present  the  results  obtained  with  the  exact  far-field  ABCs, 
both  convergence  history  and  surface  pressure  error,  for  the  flow  past  a  circular  cylinder.  In  Figures  3.2  and 
3.3  we  present  similar  results  obtained  with  the  global  and  local  far-field  ABCs,  respectively. 

From  Figures  3.1(a)  and  3.2(a)  one  can  easily  conclude  that  for  the  exact  and  global  ABCs  we  have  been 
able  to  obtain  the  multigrid  convergence  rate  according  to  the  theoretical  prediction  —  roughly  one  order  of 
magnitude  reduction  of  the  pressure  residual  per  one  W(2,l)  multigrid  cycle.  As  concerns  local  ABCs,  from 
Figure  3.3(a)  we  see  that  the  same  is  true  only  for  the  sufficiently  large  grids,  129  x  129  and  129  x  65,  and 
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consequently,  large  domains,  see  Table  3.1,  whereas  for  smaller  domains  the  convergence  slows  down. 


Circular  cylinder,  exact  boundary  condition  for  pressure 


W(2,1)  cycles 


(a)  Convergence  history 


(b)  Surface  error 


Fig.  3.1.  Computation  of  the  flow  past  a  circular  cylinder  using  exact  ABCs. 


Circular  cylinder,  global  boundary  condition  for  pressure 
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Circular  cylinder,  global  boundary  condition  for  pressure 
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(a)  Convergence  history  (b)  Surface  error 

Fig.  3.2.  Computation  of  the  flow  past  a  circular  cylinder  using  global  ABCs. 


To  assess  the  accuracy  we  compare  the  profiles  of  the  relative  error  in  pressure  (we  can  calculate  it 
explicitly  as  the  exact  solution  is  known)  on  the  surface  of  the  cylinder.  Figures  3.1(b)  and  3.2(b)  show 
that  for  the  exact  and  global  ABCs  this  relative  error  is  always  below  the  10~2  level,  i.e.,  1%,  for  all  grids. 
For  the  local  ABCs  the  accuracy  level  better  than  1%  can  be  obtained  again  only  for  the  two  largest  grid 
dimensions,  129  x  129  and  129  x  65.  With  the  reduction  of  the  grid  dimension  and,  accordingly,  domain  size 
(see  Table  3.1),  the  accuracy  deteriorates  and  eventually  reaches  the  value  0(1)  for  the  smallest  grid  129  x  9, 
see  Figure  3.3(b).  Besides,  Figure  3.3(b)  shows  that  for  the  smaller  domains  (129  x  17  and  129  x  9)  the  flow 
develops  some  spurious  non-physical  asymmetries.  Indeed,  the  two  close  lines  of  the  same  dash  pattern  on 
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Figure  3.3  (b)  correspond  to  the  error  profiles  on  the  upper  and  lower  surfaces  of  the  cylinder  for  a  particular 
grid.  Theoretically,  these  two  lines  should  be  indistinguishable,  which  is  the  case  for  the  larger  domains;  the 
actual  discrepancy  that  we  observe  for  the  smaller  domains  indicates  the  presence  of  asymmetry. 


Fig.  3.3.  Computation  of  the  flow  past  a  circular  cylinder  using  local  ABCs. 


Summarizing  for  this  case,  we  see  that  the  exact  and  global  ABCs  perform  equally  well  for  domains  of 
all  sizes,  whereas  the  performance  of  the  local  ABCs  noticeably  degrade  with  the  shrinkage  of  the  domain. 

3.3.2.  Karman-TrefFtz  airfoil.  In  Figure  3.4  we  present  computational  results  obtained  with  the 
exact  far-field  ABCs,  both  convergence  history  and  surface  pressure  error,  for  the  non-lifting  flow  past  a 
Karman-TrefFtz  airfoil  with  the  10°  trailing  edge  angle.  In  Figures  3.5  and  3.6  we  present  similar  results 
obtained  with  the  global  and  local  far-field  ABCs,  respectively. 
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Fig.  3.5.  Computation  of  the  flow  past  an  airfoil  using  global  ABCs. 


As  concerns  multigrid  convergence  rate  in  this  case,  from  Figures  3.4(a)  and  3.6(a)  we  see  that  for  the 
exact  and  local  ABCs,  respectively,  we  have  been  able  to  recover  it  according  to  the  theoretical  prediction 
(Section  3.2).  The  multigrid  convergence  rate  for  global  ABCs,  see  Figure  3.5(a),  is  also  obtained  according 
to  the  theoretical  prediction  —  about  one  order  of  magnitude  reduction  of  the  residual  level  per  one  W(2,l) 
multigrid  cycle,  but  for  smaller  grids  the  residual  never  reaches  the  machine  zero.  Instead,  starting  from  some 
level  specific  for  each  grid,  the  convergence  curve  “flattens”  and  the  residual  does  not  drop  any  further,  see 
Figure  3.5(a).  We  postpone  the  discussion  of  why,  in  our  opinion,  this  phenomenon  occurs  till  Section  3.3.3. 
Now  we  only  emphasize  that  the  level  at  which  the  residual  flattens,  see  Figure  3.5(a),  is  well  below  (several 
orders  of  magnitude)  the  truncation  error  level  for  all  the  cases  that  we  have  studied.  Therefore,  we  do  not 
expect  that  this  phenomenon  will  affect  the  final  accuracy  of  the  solution  in  any  respect. 
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Indeed,  for  the  exact  boundary  conditions  the  accuracy  for  all  grids  is  about  ~  10~4,  see  Figure  3.4(b), 
except  in  the  small  areas  near  the  leading  and  trailing  edges.  In  these  areas  we  are,  in  fact,  encountering  the 
well  recognized  problems  with  the  approximation.  Indeed,  at  the  trailing  edge,  the  conformal  mapping  that 
we  use  for  generating  the  grid  has  a  singularity,  therefore  the  surface  normal  is  not  well  defined  there.  Instead 
of  accounting  for  this  singularity  analytically  when  developing  the  discretization  we  effectively  “ignore”  it 
and  apply  the  nontranspiration  boundary  condition  by  evaluating  the  surface  normal  at  the  trailing  edge 
in  the  same  way  as  we  do  over  the  rest  of  the  of  the  airfoil  surface,  by  finite  differences  of  the  surface 
coordinates.  At  the  leading  edge,  the  numerical  viscosity  apparently  dominates  the  behavior  locally,  so  there 
is  an  artificial  dissipative  error  that  corrupts  the  total  pressure.  However,  the  latter  problems  have  no  direct 
relation  to  the  treatment  of  the  outer  boundary,  which  is  the  subject  of  the  current  study. 

Comparing  the  accuracy  provided  by  the  exact  ABCs  to  the  one  that  we  obtain  with  the  global  ABCs, 
i.e.,  comparing  Figures  3.4(b)  and  3.5(b),  we  see  that  for  the  first  three  grids,  129  x  129,  129  x  65,  and 
129  x  33,  global  ABCs  provide  for  the  same  accuracy  as  the  exact  ones  do,  i.e.,  about  10"4.  For  the  next 
smaller  grid,  129  x  17,  the  accuracy  with  the  global  ABCs  is  ~  10-3,  and  for  the  smallest  grid  129  x  9  this 
accuracy  is  ~  10“2.  Thus,  it  turns  out  that  in  terms  of  accuracy  global  ABCs  perform  somewhat  worse  than 
the  exact  ones  on  small  computational  domains  (see  Table  3.1)  for  this  airfoil  case.  It  is  clear,  however,  that 
the  accuracy  that  we  do  recover  with  the  global  ABCs,  see  Figure  3.5(b),  is  still  acceptable  for  all  purposes. 
We  also  note  that  the  exact  ABCs  are  obviously  not  available  in  realistic  situations,  while  global  ABCs  can 
be  constructed,  e.g.,  on  the  basis  of  the  difference  potentials  method,  see  [21]. 

One  might  think  that  the  slowdown  of  convergence  that  we  observe  for  global  ABCs  on  the  smaller  grids, 
see  Figure  3.5(a),  and  the  slight  deterioration  of  accuracy  that  we  also  observe  for  these  boundary  conditions 
on  the  smaller  grids,  see  Figure  3.5(b),  could  be  related  to  one  another.  We  claim  that  this  is,  in  fact,  not 
so,  because  for  the  local  ABCs  that  do  provide  in  this  case  for  the  optimal  convergence  till  the  machine  zero 
on  all  grids,  see  Figure  3.6(a),  the  accuracy  rapidly  deteriorates  as  the  domain  shrinks,  see  Figure  3.6(b), 
and  reaches  the  level  ~  10-1,  i.e.,  10%,  for  the  smallest  grid  129  x  9. 

Summarizing  for  the  airfoil  case,  we  see  that  the  global  ABCs  actually  perform  only  slightly  worse  than 
the  exact  ones  and  their  accuracy  is  certainly  within  the  practically  acceptable  limits  on  all  grids,  whereas 
the  solutions  obtained  with  the  local  ABCs  rapidly  lose  accuracy  as  the  domain  size  reduces. 

3.3.3.  Comparison  of  the  cylinder  and  airfoil  computations.  The  first  easy  observation  that  one 
can  make  is  that  on  the  average  the  relative  accuracy  for  the  Karman-Trefftz  airfoil  is  better  than  that  for 
the  circular  cylinder  for  all  the  computations  that  we  have  conducted,  see  Figures  3.1(b)  —  3.6(b).  With  no 
rigorous  explanation  we  can  attribute  this  to  the  mere  fact  that  in  some  sense  a  “thick”  cylinder  introduces 
more  of  a  perturbation  into  the  flow  that  a  “thin”  airfoil  does  (the  thickness  ratio  of  the  airfoil  is  15%, 
see  Section  3.1).  The  same  fact  apparently  accounts  for  the  slowdown  of  multigrid  convergence  on  smaller 
domains  when  computing  the  flow  past  a  cylinder  with  local  ABCs,  see  Figure  3.3(a),  whereas  the  same 
local  ABCs  for  the  airfoil  still  provide  for  the  optimal  multigrid  convergence  rate,  see  Figure  3.6(a).  We 
emphasize  here  that  the  convergence  slowdown  on  Figure  3.3(a)  (local  ABCs  for  the  cylinder)  is  different  in 
nature  to  the  one  on  Figure  3.5(a)  (global  ABCs  for  the  airfoil)  because  in  the  latter  case  we  do  recover  the 
optimal  convergence  rate  and  lose  it  only  toward  the  end  of  the  computation,  when  the  residuals  are  already 
sufficiently  small,  and  in  the  former  case  the  convergence  is  suboptimal  from  the  very  beginning. 

The  phenomenon  of  “flattening”  of  the  residuals  for  global  ABCs,  see  Figure  3.5(a),  can  apparently 
be  attributed  to  the  following.  The  boundary  conditions  have  actually  been  constructed  on  the  conformal 
plane  £,  where  the  body  is  always  a  cylinder  and  the  grid  is  always  exactly  polar,  which,  in  particular, 
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means  orthogonal.  In  Sections  2.1  and  2.2.1  we  have  addressed  the  importance  of  maintaining  the  exact 
conservation  on  the  grid  from  the  standpoint  of  solvability  and  also  shown  (see  equations  (2.40)  and  (2.43)) 
how  to  maintain  the  exact  conservation  on  this  particular  grid,  i.e.,  grid  of  polar  coordinates.  For  the  flow 
around  the  airfoil  the  discretization  is  built  directly  using  finite  volumes  on  the  actual  grid  in  the  physical 
plane.  For  small  domains,  this  grid  may  have  slight  numerical  deviations  from  orthogonality  even  near  the 
outer  boundary.  This,  in  turn,  results  in  a  slightly  different  discrete  operator.  (As  has  been  mentioned  in 
Section  2.2.1,  in  contrast  to  (2.28)  it  may  contain  diagonal  nodes  with  small  coefficients.)  This  may  give  rise 
to  a  minor  violation  of  the  overall  conservation  on  the  grid,  the  smaller  the  domain,  the  larger  the  violation. 
(Indeed,  the  left-hand  side  of  (2.40)  and  (2.43)  has  to  have  this  particular  form  to  be  incorporated  into 
the  boundary  conditions,  and  the  conservation  will  now  be  achieved  with  a  different  operator.)  This  small 
“incompatibility”  is  responsible  for  the  flattening  of  the  convergence  history  curve  (i.e.,  residual  level  vs. 
number  of  multigrid  cycles,  see  Figure  3.5(a)).  A  similar  phenomenon  was  observed  in  [13]  when  solving  a 
Neumann  problem,  for  which  the  solvability  condition  was  satisfied  exactly  for  the  continuous  formulation, 
while  only  up  to  the  truncation  error  level  for  the  discrete  formulation. 

3.4.  Solving  the  Euler  system.  Let  us  first  note  that  in  the  framework  of  the  full  Euler  system  we 
are  computing  only  the  flow  around  the  airfoil.  The  cylinder  case  analyzed  in  Section  3.3  can,  in  fact,  be 
used  only  for  the  pure  pressure  computations.  It  is  not  suitable  for  the  full  Euler  equations  because  the 
artificial  dissipation  of  the  scheme  will  result  in  flow  separation.  This  effect  is,  again,  due  to  the  “thickness” 
of  the  cylinder. 

In  Figure  3.7  we  present  convergence  histories  on  all  grids  for  all  three  types  of  the  ABCs.  As  one  can 
see,  the  convergence  in  all  cases  is  somewhat  slower  than  the  theoretically  predicted  rate  of  0.125  per  one 
W(2,l)  cycle  that  we  use.  The  slowdown  of  the  multigrid  convergence  that  we  observe  here  is  obviously 
not  related  to  the  treatment  of  the  artificial  external  boundary  because  it  takes  place  for  all  three  types 
of  the  ABCs.  As  mentioned  in  Section  3.2,  this  slowdown  is  most  likely  related  to  the  apparently  invalid 
omission  of  subprincipal  terms  near  the  stagnation  point  in  the  construction  of  the  relaxation  scheme.  This 
issue  is  discussed  in  [16].  We  also  emphasize  that  for  the  local  ABCs,  see  Figure  3.7(c),  the  convergence 
is  noticeably  slower  than  for  the  other  two  types,  see  Figures  3.7(a)  and  3.7(b),  and  for  the  smallest  grid 
129  x  9  the  solution  with  local  ABCs  does  not  converge  at  all.  That’s  why  we  have  only  four  curves  (as 
opposed  to  five)  plotted  on  Figure  3.7(c). 

In  Figure  3.8  and  3.9,  we  are  showing  the  surface  pressure  error  profiles  for  all  five  grids  and  three  types 
of  the  ABCs  that  we  have  used.  As  has  been  mentioned,  the  solution  with  local  ABCs  on  the  smallest  grid 
129  x  9  did  not  converge,  therefore  the  corresponding  error  curve  is  not  available.  (Note,  as  opposed  to 
Section  3.3,  here  we  present  separate  error  profile  plots  for  each  grid  and  all  three  types  of  the  ABCs,  i.e., 
have  three  curves  per  plot,  rather  than  showing  five  curves  from  all  five  grids  on  one  plot  for  each  geometry 
and  each  type  of  ABCs.) 

Figures  3.8  and  3.9  allow  one  to  conclude  that  for  the  exact  ABCs  the  surface  pressure  error  is  always 
below  the  10-2  level,  for  global  ABCs  it  is  also  below  than  10~2,  except  for  the  smallest  domain,  for  which 
it  is  only  slightly  above  the  level  of  10-2,  i.e.,  1%.  For  local  ABCs  the  accuracy  rapidly  deteriorates  with 
the  domain  size  decrease  and  eventually  reaches  0(1)  on  the  grid  of  moderate  dimension  129  x  17  before  the 
solution  breaks  on  the  smallest  grid  129  x  9. 

As  opposed  to  Section  3.3,  here  we  are  computing  all  three  flow  quantities  rather  than  only  pressure. 
Therefore,  we  also  compare  error  profiles  for  velocities.  In  Figure  3.10  we  present  velocity  errors  for  the 
smallest  grid,  on  which  all  three  solutions  converged,  i.e.,  for  the  grid  129  x  17.  Figure  3.10(a)  shows  that 
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Fig.  3.7. 


Convergence  histories  for  the  flow  past  an  airfoil  computed  in  the  framework  of  the  full  Euler  equations . 


both  the  exact  and  global  ABCs  perform  equally  well  for  the  surface  velocity  (the  error  is  slightly  above 
10“2)  and  the  solution  obtained  with  local  ABCs  substantially  lacks  accuracy.  At  the  outer  boundary,  see 
Figure  3.10(b),  global  ABCs  provide  for  the  error  level  below  10-2,  whereas  local  ABCs  are  lagging  behind 
by  roughly  two  orders  of  magnitude. 

Note,  we  do  not  present  the  error  curve  for  the  exact  ABCs  on  Figure  3.10(b)  because  the  Dirichlet  data 
for  these  boundary  conditions  are  taken  from  the  exact  solution  and  thus  the  corresponding  error  would  have 
been  identically  equal  to  zero. 

Summarizing  for  the  full  Euler  airfoil  case,  we  see  that  global  ABCs  are  almost  as  good  as  the  exact  ones 
in  terms  of  both  the  solution  accuracy  and  multigrid  convergence  rate.  As  concerns  global  vs.  local  ABCs, 
from  Figures  3.7,  3.8,  3.9,  and  3.10  we  conclude  that  global  ABCs  clearly  outperform  the  local  ones  from 
the  standpoints  of  accuracy,  convergence  rate,  and  robustness.  Accuracy- wise,  they  allow  to  reduce  the  grid 
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dimension  by  a  factor  of  8-16  with  no  or  very  little  increase  in  error  and  with  the  obvious  corresponding 
reduction  in  computational  costs  —  this  is  the  primary  result  that  we  have  expected.  Regarding  the  robust¬ 
ness,  global  ABCs  provided  for  the  convergence  on  the  smallest  grid  129  x  9,  when  the  algorithm  with  local 
ABCs  simply  failed  to  converge. 


Surface  pressure  error  for  Karman-Trefftz  airfoil, 
10°  trailing  edge  angle,  129  x  129  grid 


Surface  pressure  error  for  Karman-Trefftz  airfoil, 
10°  trailing  edge  angle,  129  x  65  grid 


Surface  pressure  error  for  Karman-Trefftz  airfoil, 
10°  trailing  edge  angle,  129  x  33  grid 


Surface  pressure  error  for  Karman-Trefftz  airfoil, 
1 0°  trailing  edge  angle,  1 29  x  1 7  grid 


Fig.  3.8.  Surface  pressure  error  profiles  for  the  flow  past  an  airfoil  computed  in  the  framework  of  the  full  Euler  equations. 


4.  Discussion  and  conclusions.  Motivated  by  the  two  independent  successful  developments  in  CFD 
that  have  appeared  recently:  The  new  factorizable  schemes  for  the  equations  of  hydrodynamics  that  facilitate 
the  construction  of  optimally  convergent  multigrid  algorithms,  and  highly  accurate  global  far-field  artificial 
boundary  conditions,  we  have  built  and  tested  a  unified  methodology  for  calculating  incompressible  flow 
solutions  based  on  the  combination  of  the  two  aforementioned  approaches.  The  primary  result  that  we  have 
obtained  is  the  following.  Global  ABCs  do  not  hamper  the  optimal  (i.e.,  unimprovable)  multigrid  convergence 
rate  pertinent  to  the  solver.  At  the  same  time,  contrary  to  the  standard  local  ABCs,  the  solution  accuracy 
provided  by  the  global  ABCs  deteriorates  very  slightly  or  does  not  deteriorate  at  all  when  the  computational 
domain  shrinks,  which  clearly  translates  into  substantial  savings  of  computer  resources. 
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Fig.  3.9.  Surface  pressure  error  profiles  for  the  flow  past  an  airfoil  computed  in  the  framework  of  the  full  Euler  equations . 


Speed  error  at  surface  for  Karma n-T re fftz  airfoil, 
10°  trailing  edge  angle,  129  x  17  grid 
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Global  boundary  condition 
Local  boundary  condition 
Exact  DIrichlet  boundary  condition 


Speed  error  at  outer  boundary  for  Karman-Trefftz  airfoil, 
1 0°  trailing  edge  angle,  129x17  grid 


(a)  On  the  surface 


(b)  On  the  outer  boundary 


Fig.  3.10.  Velocity  error  profiles  for  the  flow  past  an  airfoil  computed  in  the  framework  of  the  full  Euler  equations  on  the 
129  x  17  grid. 


The  combined  methodology  was  developed  for  the  most  simple  formulation  of  both  the  factorizable 
scheme  and  the  ABCs:  the  so-called  pressure  Poisson  scheme,  and  the  ABCs  built  semi-analytically  using 
the  separation  of  variables  along  the  artificial  boundary.  The  methodology  was  tested  on  a  class  of  incom¬ 
pressible  inviscid  non-lifting  two-dimensional  flows  with  the  exact  solutions  readily  available  through  the 
implementation  of  the  classical  complex  variable  technique.  On  one  hand,  this  intentionally  elementary  for¬ 
mulation  (both  theoretically  and  experimentally)  allows  to  use  the  analytical  approaches  and  results  as  much 
as  possible  for  both  constructing  the  methodology  and  numerically  assessing  its  performance.  On  the  other 
hand,  this  formulation  still  keeps  the  key  fundamental  principles,  on  which  both  the  factorizable  scheme 
and  global  ABCs  are  based,  and  thus  should  be  regarded  not  as  an  isolated  model  example  but  rather  as  a 
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foundation  of  the  future  approaches  to  more  complex  problems. 

In  a  series  of  numerical  tests  (Section  3)  we  have  compared  three  types  of  the  ABCs:  Exact  Dirichlet 
ABCs  obtained  from  the  exact  analytical  solution,  global  ABCs,  and  standard  local  ABCs  of  the  kind  that 
is  used  extensively  and  routinely  nowadays  when  neither  the  exact  data  nor  any  advanced  methodology  is 
available.  From  our  computations  we  can  conclude  that  convergence- wise,  global  ABCs  perform  practically 
as  good  as  the  exact  ones;  in  other  words  they  recover  the  optimal  multigrid  convergence  rate  that  cannot  be 
improved.  From  the  standpoint  of  accuracy,  they  allow  for  a  very  substantial  reduction  of  the  computational 
domain  size,  two  orders  of  magnitude  in  terms  of  the  actual  diameter  and  at  least  one  order  of  magnitude  in 
terms  of  the  grid  dimension,  with  no  or  very  little  deterioration  in  the  final  quality  of  the  calculated  solution. 
This  obviously  yields  a  considerable  reduction  in  the  corresponding  computational  costs,  which  altogether 
constitutes  the  primary  result  that  we  have  expected. 

The  aforementioned  reduction  in  the  computational  domain  size  is  obtained  when  comparing  the  per¬ 
formance  of  the  global  ABCs  with  that  of  the  standard  local  boundary  conditions.  We  should  emphasize 
here  that  as  the  exact  ABCs  are  not  attainable  routinely,  the  comparative  assessment  of  global  and  local 
methodologies  is  of  the  foremost  importance  from  the  viewpoint  of  computational  practice.  Our  numerical 
experiments  show  that  besides  the  solution  accuracy,  global  ABCs  clearly  surpass  the  local  ones  from  the 
standpoints  of  multigrid  convergence  rate  and  robustness.  (The  latter  point  refers  to  the  case  when  global 
ABCs  could  provide  for  convergence  while  the  local  ones  diverged.)  A  similar  kind  of  behavior  has  been 
observed  previously  for  the  global  DPM-based  ABCs  combined  with  the  older  suboptimal  multigrid  method¬ 
ologies  [21,25-27,30-33].  As  at  the  moment  different  flavors  of  local  boundary  conditions  still  dominate  the 
area  of  production  computations  in  CFD,  the  new  combined  technique  developed  and  tested  in  this  paper 
provides  a  potential  for  creating  new  flow  solvers  far  superior  to  those  currently  in  use. 

An  important  implementation  lesson  learned  when  building  the  unified  methodology  is  that  the  inho¬ 
mogeneities  in  the  elliptic  factor  of  the  system  have  to  be  handled  very  carefully.  In  particular,  the  exact 
conservation  on  the  grid  has  to  be  strictly  enforced  to  guarantee  the  solvability  of  the  discrete  problem. 
Let  us  note  that  the  actual  inhomogeneity  of  the  pressure  Poisson  equation  that  we  have  studied  in  this 
paper  is  not  going  to  be  encountered  in  the  future  because  the  forthcoming  fact oriz able  scheme  for  the 
incompressible  case  will  be  based  on  the  Laplace  equation  for  the  velocity  potential  as  an  elliptic  factor 
(always  homogeneous).  On  the  other  hand,  the  current  study  provides  foundation  for  treating  the  compress¬ 
ible  flow  equations,  in  which  the  full-potential  part  will  always  involve  far-field  inhomogeneities  due  to  the 
sub- principal  terms.  The  boundary  conditions  for  this  equation  will  be  the  same  as  those  we  have  studied 
here:  Neumann’s  boundary  condition  on  the  solid  wall  and  requirement  of  boundedness  of  the  solution  at 
infinity.  Therefore,  the  conclusions  of  the  current  study,  in  particular  those  related  to  the  solvability  issues, 
are  going  to  be  useful. 

Our  future  plans  involve  extensions  along  both  experimental  and  theoretical  lines.  Numerically,  we  will 
incorporate  a  different  formulation  of  the  scheme  and  study  more  complex  cases,  including  compressible 
flows.  Theoretically,  we  will  try  and  analyze  the  properties  of  the  discrete  ABCs’  operator  (a  generalization 
of  the  operator  T,  see  (2.35)  and  (2.37))  in  terms  of  the  influence  it  may  exert  on  multigrid  convergence, 
e.g.,  whether  the  matrix  of  the  corresponding  overall  system  that  is  solved  by  iterations  may  appear  either 
symmetric  positive  definite  or  an  M-matrix  (see  [41]  for  the  definition  of  the  latter).  Global  ABCs  for 
the  more  complex  cases  that  we  plan  on  studying  in  the  future  will  be  obtained  using  difference  potentials 
method. 
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